yProcessingClub

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

調布祭に行くなどした

2017年11月25日、母校の電気通信大学の学園祭「調布祭」に行ってきた。

 

久々に調布を訪れたのだが、駅前が栄えていてたまげた。

写真には撮っていないが、ビックカメラがオープンしていた。

 

 

電通大はお変わりないようで。

 

 

テキトーに屋台のクレープなんかを食べる。

 

目的は学園祭そのものというより、大学時代の友達と会うことにあった。

友達と合流したので大学を早々に後にし、ご飯を食いに行く。

3人とも1ポンドステーキを頼む。味はまぁまぁという感じだった。 

 

 

ご飯屋さんを探していて、ラーメン屋「そらまめ」が閉店していたことを知って若干悲しむ。

「そらまめ」は知り合いの間では好評であった。(自分としてはそこまで美味しいとは思わなかったが・・・)

 

 

 

久々に会った友達はこの後サークル?のOB会に参加するというので離脱し、残った二人で秋葉に行った。

 

この友達とはよく秋葉に遊びに行くが、毎回特に何をするということも無くただステーキを食って終わりなのである。

この日はすでに調布でステーキを食べたので本当に何もせず解散した。何で秋葉に行った?

 

 

 

池袋のジュンク堂へ。

 

もともとジュンク堂に行く予定だったのだが、遠方に就職してなかなか会えない友達が調布祭に来てたというので急遽予定を変更したのである。一人になったので本来の予定を回収していく。

 

 

 

 信号処理の超基礎的な本とベイズ統計の本を買った。

せっかく来たんだしもっと色々回りたかったが、体力が限界で目眩がしてきたので退散。

 

 

 ジュンク堂の近くのスシローで晩飯。

 

 

疲れたけど楽しかった(こなみ)

 

味奈登庵 富士山盛り

味奈登庵 富士山盛り。 麺は1kgだそうである。

天ぷらもつけたので相当なボリュームだと覚悟していたが、 神豚の大豚ダブルなど日常的に食べている身からすると余裕であった。

蕎麦や天ぷら自体のクオリティも悪くなく、普通に美味しかった。

http://www.minatoan.com/menu/menu_ss/ss_main/ss_main.html

肉饂飩とみ坂 ピリ辛冷やし肉うどん

うまかったー

量がまったくもって足りないのでお店のハシゴ不可避であった。

食べログはここ↓
https://tabelog.com/kanagawa/A1401/A140101/14062998/

matlab スペクトログラム 周波数軸および時間軸方向に切り取る

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の値を使用すると良いかと思われる.
後ほど修正予定.
申し訳ないが取り急ぎ報告のみ.


スペクトログラムは横軸時間,縦軸周波数,レベルを色の濃淡で表したグラフである.

任意の周波数で横方向に切り取ったグラフおよび任意の時間で縦方向に切り取ったグラフを描画してみる.


図解するとこんな感じ.

今回はこれらを関数として用意し,mainで呼び出すことにした.

main

% 分析信号としてチャープ信号に白色ノイズを重畳したものを用いる

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

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

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


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

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

% STFTする
Nw = 512;                   % 時間窓の長さ
Nol = round(Nw/2);          % オーバーラップ
NFFT = Nw;                  % FFT点数
[~, freq, time, ps] = spectrogram(z,Nw,Nol,NFFT,Fs,'yaxis');


% 周波数-レベルグラフを描く
t = 2;
figure;
plot_freq_level_plane_of_spectrogram(ps, Fs, t);

% 時間-レベルグラフを描く
f = 10*10^3;
figure;
plot_time_level_plane_of_spectrogram(ps, Fs, f);

plot_freq_level_plane_of_spectrogram()

function [] = plot_freq_level_plane_of_spectrogram(ps,Fs,t)
% スペクトログラムの時間t(秒)における周波数-レベル平面を描画する

% 周波数配列の計算
N = size(ps,1);                 % 周波数方向のサイズ
freq = ((1:N)-1)*Fs/N/2/1000;   % 周波数配列(kHz単位)

% 指定時間(秒)における時間-レベルの配列を取得
d_freq = Fs / N;                % 周波数分解能
d_time = 1 / d_freq;            % 時間分解能
time = round(t / d_time);       % 時間軸インデックスを計算
P = ps(:,time);                 % 配列のコピー

% グラフを描画する
P = 10*log10(P);                % 振幅のdBへの変換
plot(freq,P);

% グラフの描画設定
xlabel('周波数(kHz)');        % X軸ラベル
xlim([0, Fs/2/1000]);           % X軸範囲(ナイキスト周波数(kHz)まで)
ylabel('パワー/周波数(dB/Hz)')  % Y軸ラベル
ylim([-60, 0]);                 % Y軸範囲
str = sprintf('周波数-レベルグラフ(t=%0.1f秒)',t);
title(str);                     % タイトル

end

plot_time_level_plane_of_spectrogram()

function [] = plot_time_level_plane_of_spectrogram(ps,Fs,f)
% スペクトログラムの周波数f(Hz)における時間-レベル平面を描画する

% 時間軸配列を計算
N = size(ps,1);                 % 周波数軸方向のサイズ
d_freq = Fs / N;                % 周波数分解能
d_time = 1 / d_freq;            % 時間分解能
N = size(ps,2);                 % 時間軸方向のサイズ
time = ((1:N)-1) * d_time;      % 時間軸配列

% 指定周波数(Hz)における時間-レベルの配列を取得
freq = round(f / d_freq * 2);   % 周波数インデックスを計算
P = ps(freq,:);                 % 配列をコピーする

% グラフを描画する
P = 10*log10(P);                % 振幅のdBへの変換
plot(time,P);

% グラフの描画設定
xlabel('時間(secs)');            % X軸ラベル
xlim([0, (N-1)*d_time]);         % X軸範囲
ylabel('パワー/周波数(dB/Hz)')   % Y軸ラベル
ylim([-60, 0]);                  % Y軸範囲
str = sprintf('時間-レベルグラフ(f=%0.1fkHz)',f/1000);
title(str);                      % タイトル

end

結果


fig.周波数-レベルグラフ(t = 2秒)


fig.時間-レベルグラフ(f = 10kHz)



期待する位置にピークが出ており,正しくグラフが描けたと思われる.


 

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

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で描画したグラフ(いろいろ設定後)



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



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