yProcessingClub

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

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

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

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

f:id:Yuri-Processing-Club:20170916102546j:plain
図解するとこんな感じ.

今回はこれらを関数として用意し,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

結果

f:id:Yuri-Processing-Club:20170916102920p:plain
fig.周波数-レベルグラフ(t = 2秒)

f:id:Yuri-Processing-Club:20170916100604p:plain
fig.時間-レベルグラフ(f = 10kHz)



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