yProcessingClub

すみません、許してください

matlab スペクトログラム をmeshで描画してみる

matlab

jp.mathworks.com

matlabにはspectrogramという関数があり,これを使うとスペクトログラムが描ける.
今回はmesh関数で同じように描けるか試してみる.

2022/03/29 追記：重大なバグを発見した．
[~, freq, time, ps] = spectrogram(z,Nw,Nol,NFFT,Fs,'yaxis');
とあるが，psは確率密度関数なのでspectrogramの結果を使用していない．
[s, freq, time, ~] = spectrogram(z,Nw,Nol,NFFT,Fs,'yaxis');
としてpsのかわりにsの値を使用すると良いかと思われる．
文中で「色合いが若干違うのはなんだろ・・・計算等や使っている変数が違うかもしれないので,要調査.」と書いているが，この辺のバグのせいかと思われる．
後ほど修正予定．
申し訳ないが取り急ぎ報告のみ．

準備

% 分析する信号zを生成する

Fs = 44100;                 % サンプリング周波数(Hz)
T = 10;                     % 信号の時間長さ（秒）
N = Fs * T;                 % 信号のサンプル数

n = randn(1, N);            % N点の白色ノイズ

% 線形チャープ信号の生成
f1 = 5000;                  % 開始周波数(Hz)
f2 = 15000;                 % 終了周波数(Hz)
t = ((1:N) - 1) / Fs;       % 周期1/Fs, サンプル数Nの時間配列
s = chirp(t,f1,T,f2);       % チャープ信号の生成

% ノイズとチャープ信号を足し合わせる
SNRdB_original = snr(s,n);  % 元のSN比を求める
SNRdB = -10;                % 設定するSN比
coef = 10^((SNRdB_original - SNRdB)/20);    % SN比を調節するための係数
n = coef * n;               % ノイズパワーを調整

z = s + n;                  % チャープ信号とノイズを重畳する

分析する信号はチャープ信号に白色ノイズを重畳したものにした.



spectrogramで描画

% STFTする

Nw = 512;                   % 時間窓の長さ
Nol = round(Nw/2);          % オーバーラップ
NFFT = Nw;                  % FFT点数
spectrogram(z,Nw,Nol,NFFT,Fs,'yaxis');  % スペクトログラムを描画
caxis([-60 0]);             % レベルの範囲を指定
colormap jet;               % カラーマップを設定


fig.spectrogramで描画したグラフ

これをmeshで再現してみる.

jp.mathworks.com



meshを使って描画(修正前)

% meshを使って描画する

[~, freq, time, ps] = spectrogram(z,Nw,Nol,NFFT,Fs,'yaxis');
psdB = 10*log10(ps);        % dBに変換
figure;                     % figureウィンドウの作成
mesh(time, freq, psdB);     % パワースペクトルを描画
caxis([-60 0]);             % レベルの範囲を指定
colormap jet;               % カラーマップを設定


fig.meshをそのまま使ったグラフ

視点や軸ラベルが設定されていないので,これを設定する.


meshにいろいろ設定する

% meshを使って描画する（修正後）

[~, freq, time, ps] = spectrogram(z,Nw,Nol,NFFT,Fs,'yaxis');
psdB = 10*log10(ps);        % dBに変換
figure;                     % figureウィンドウの作成
freq = freq / 1000;         % HzをkHzに変換
mesh(time, freq, psdB);     % パワースペクトルを描画
caxis([-60 0]);             % レベルの範囲を指定
colormap jet;               % カラーマップを設定
view(0,90);                 % 見る角度を指定
xlabel('時間(secs)');       % X軸ラベル
ylim([0,Fs/2/1000]);        % Y軸範囲を指定(ナイキスト周波数(kHz)まで)
ylabel('周波数(kHz)');      % Y軸ラベル
h = colorbar;               % カラーバーを表示
ylabel(h, 'パワー/周波数(dB/Hz)');% カラーバーのラベル


fig.meshで描画したグラフ(いろいろ設定後)


結果

上で貼ったグラフを再度貼り付ける.


fig.spectrogramで描画したグラフ


fig.meshで描画したグラフ(いろいろ設定後)



ほぼ同一のグラフを描くことができた.



色合いが若干違うのはなんだろ・・・
計算等や使っている変数が違うかもしれないので,要調査.


　

matlab SN比計算　SN比から信号の大きさを逆算する

matlab

yuri-processing-club.hatenablog.com

こちらの記事でSN比を指定してノイズの大きさを調整する計算を紹介した。
matlabにて計算が正しいかを確認する。

matlabコード

% sin信号及び白色ノイズを生成
% 自前の計算及び標準の関数を用いてSN比を計算し、表示する
% 指定したSN比になるようノイズのパワーを調整する

Fs = 44100;                 % サンプリング周波数(Hz)
T = 10;                     % 信号の時間長さ（秒）
N = Fs * T;                 % 信号のサンプル数

n = randn(1, N);            % N点の白色ノイズ

% sin信号の生成
f = 100;                    % sin波の周波数(Hz)
t = ((1:N) - 1) / Fs;       % 周期1/Fs, サンプル数Nの時間配列
s = sin(2 * pi * f * t);    % N点, 周波数fのsin波

% パワーの計算
Pn = n * n';                % ノイズのパワー(2乗の和)
Ps = s * s';                % チャープ信号のパワー(2乗の和)

% SN比の計算
SNR = 10 * log10(Ps/Pn);
SNR1 = snr(s, n);           % matlab標準の関数
fprintf('調整前：SNR = %fdB, SNR1 = %fdB\n', SNR, SNR1);

% 指定したSN比となるようにノイズの大きさを調整する
SNR2 = -30;                 % 指定するSN比(dB)
coef = 10^((SNR - SNR2)/20);% ノイズの大きさを調整する係数
n = coef * n;               % ノイズの大きさを調整

% 再度SN比の計算
Pn = n * n';                % ノイズのパワーを再度求める
SNR = 10 * log10(Ps/Pn);
SNR1 = snr(s, n);
fprintf('調整後：SNR = %fdB, SNR1 = %fdB\n', SNR, SNR1);

実行結果

調整前：SNR = -3.006121dB, SNR1 = -3.006121dB
調整後：SNR = -30.000000dB, SNR1 = -30.000000dB


となり、正しい結果となった。




　

matlabを導入

matlab

仕事の研究でmatlabを使っている。
家でも色々やりたいと思って自費で購入することにした。

f:id:Yuri-Processing-Club:20170908062535p:plain

バージョンは2017a。
仕事で使っているのは2014？だったと思うので、だいぶ進んでいる。


これで家でも研究が捗るというわけである。






　 　

SN比計算　SN比から信号の大きさを逆算する

信号処理

シミュレーション等で信号にノイズを重畳する際、SN比を指定したいことがあると思われる。
検索しても出てこなかったので、自分で記事を書いてみようと思う。
合ってるかは知りません。


[2018/11/11追記ここから]
合っているかは知らないと書いたが、おそらく合っていない（途中計算が）。
yuri-processing-club.hatenablog.com
ここに書いた通りであるが、離散信号の積分計算をする場合、T_sを掛ける必要がある。 ただしP_s/P_nを計算する際に分母分子にあるT_sが約分されて消えるので、T_sがあってもなくても結果は変わらないはず。 そして今回はT_s=1として正規化していることにして、以降の説明は特に修正しないものとする。
[2018/11/11追記ここまで]

SN比が任意の値になるように雑音をa倍する。
この係数aを求めるものとする。

まず、元のSN比を求めてみる。

サンプル数K個の信号s\left[k\right]及びノイズn\left[k\right]として(k=0,1,2,\cdots,K-1)、
\bf s \rm = [s [ 0 ] \ s [ 1 ] \ s [ 1 ] \ \cdots \ s [ K-1 ]]
\bf n \rm = [n [ 0 ] \ n [ 1 ] \ n [ 1 ] \ \cdots \ n [ K-1 ]]

2乗して、
\bf s^2 \rm = [s [ 0 ]^2 \ s [ 1 ]^2 \ s [ 1 ]^2 \ \cdots \ s [ K-1 ]^2]
\bf n^2 \rm = [n [ 0 ]^2 \ n [ 1 ]^2 \ n [ 1 ]^2 \ \cdots \ n [ K-1 ]^2]

各要素を足し合わせるとパワーP_s及びP_nが求まる。
P_s = \sum \bf s^2 \ \rm = \displaystyle \sum _{k=0} ^{K-1} s [ k ] ^2 = s [0 ]^2 + s [1 ]^2 + s [2 ]^2 + \cdots + s [K-1 ]^2

P_n = \sum \bf n^2 \ \rm = \displaystyle \sum _{k=0} ^{K-1} n [ k ] ^2 = n [0 ]^2 + n [1 ]^2 + n [2 ]^2 + \cdots + n [K-1 ]^2

SN比は以下の通り。

SNR_{original} = 10\log (\frac{P_s}{P_n})

次にノイズ\bf nをa倍したものを\bf n_{sc}として、\bf sと\bf n_{sc}のSN比が指定した値となるような係数aを求める。
\bf n_{sc} \rm = a \bf n

パワーP_{n_{sc}}は

P_{n_{sc}} = \sum \bf n_{sc} ^2 = \rm \sum a^2 \bf n^2 \ \rm = a ^2 \sum \bf n^2 \rm = \displaystyle a ^2 P_n


SN比は以下の通り。

SNR = 10\log(\frac{P_s}{P_{n_{sc}}}) = 10\log (\frac{P_s}{a ^2 P_n}) = 10\log (\frac{P_s}{P_n}) - 10\log a^2 \\ \ \ \ \ \ \ \ \ \ = SNR_{original} - 10\log a^2

*1

ここからaを求めていく。

10\log a^2 = SNR_{original} - SNR \\ \ \ 20\log a = SNR_{original} - SNR \\ \ \ \ \ \ \ \log a = \frac{SNR_{original} - SNR}{20}\\ \\ \ \ \ \ \ \ \ \ \ \ \ \ a = 10^{\frac{SNR_{original} - SNR}{20}}

ということで、SNRに任意の値を入れることで、係数aが求められる。

*1:はてなブログはtexのイコール揃えが出来ないっぽい？ので、スペース挿入の力技でイコールを揃えている。

がっこうぐらし考察　めぐねえの死体はどこに？

がっこうぐらし！

がっこうぐらし！考察。
（気づきの掲示がメインで、考察についてはガバガバです。）

f:id:Yuri-Processing-Club:20170820212802j:plain

地下を彷徨っていためぐねえはみーくんによって解放されるわけですが、
その後物資を求めて皆で行った時にはめぐねえの死体が消えています。

学園生活部の部室から物資があるエリアに至るまでの動線の外、 例えば地下の奥などで処理したなら、物資を取りに来た生活部メンバーが死体を見ずに済みます。
しかし、みーくんはご丁寧に地下に行く際に絶対通る場所であるシャッター前で処理しています。

だのにシャッター付近にはめぐねえの死体どころか血溜まりすらありません。

考察１　由紀の視点である

物資を取りに行った際にシャッターで見た光景は由紀の視点であり、
由紀からは何も見えていなかった。
・・・さすがの由紀といえどもめぐねえの格好をした死体が転がってたら何らかの反応をするのでは？

考察２　みーくんはめぐねえを倒していない

みーくんがシャッター前でめぐねえを倒していないので死体も無いという理屈。

みーくんの目的は薬の入手であり、めぐねえの討伐ではありません。
めぐねえから逃げながら、もしくはそもそも遭遇せずに薬を入手してもいいわけです。

りーさん「そう　やっぱり･･･　めぐねえだったんだ」
みーくん「･･････　はい･･････」

と会話してるので、地下にめぐねえがいたことは伝えています。
ただし倒したとは伝えていません。

めぐねえが倒されたシーンは読者しか見ていないので、 このシーンが嘘であったとしても、 みーくんはりーさん達には嘘をついていないことになります。

ただしその場合、物資を取りに行った際にゾンビめぐねえに襲われる心配はしなくていいのかということになります。

考察３　みーくんはシャッターの前ではない場所でめぐねえを倒した

これならばめぐねえの死体がシャッター前に転がっていない件も、 みーくんがめぐねえを倒していないとしたらゾンビめぐねえに襲われる心配をせねばいけない件にも説明がつきます。
嘘をつかれたのは読者だけですので。

考察４　めぐねえは地下にはいなかった

地下は平和そのもので、誰もいなかった。

まずくるみちゃんが遭遇しためぐねえは髪がロングです。
由紀が見ている幻覚めぐねえもそうですが、ロングめぐねえ＝幻覚の可能性があります。

f:id:Yuri-Processing-Club:20170820222835j:plain
この通り、地下に行った時点でショートですからね。

くるみちゃんは地下へ行く前、めぐねえが書類を隠し持っていたことを知ってかなり動揺していました。
よって、めぐねえの幻覚を見たとしても不思議ではありません。

f:id:Yuri-Processing-Club:20170820225529j:plain
くるみちゃんの傷跡は噛まれた痕っぽいですが、めぐねえは爪で切り裂いてます。
実際の傷と、くるみちゃんが襲われた光景に不整合があるわけです。
爪で服を切り裂いた後でガブっといったのかもしれませんが、噛まれた描写はされていません。

というわけで幻覚の可能性があります。
（じゃあくるみちゃんは一体誰に襲われたんだ？となりますが。きっとその辺のモブゾンビに襲われたのでしょう。）

次にみーくんですが、遭遇したのはショートめぐねえであり、幻覚ではないと思われます。
みーくんがショートめぐねえと遭遇したのが真実だった場合、
考察３で書いた通り、シャッター前でない場所でめぐねえを倒したか、
もしくは「めぐねえと会った」と故意に嘘をついたと考えるのが妥当ではないでしょうか。

考察５　別の時空のものが混じっている

同じ時空の場合、みーくんがりーさんに嘘をつくなどしないと成立しないように思われます。
やはり別の時空の話が混じっていると考えるべき？