yProcessingClub

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

matlab 信号の両端を滑らかに0にする場合はコサインテーパーを使う

 実験用に正弦信号を作成することがある.例として1~1.5秒部分に5kHzの正弦波,1.5~3秒部分に10kHzの正弦波が含まれる信号の時間-振幅グラフを図1に,スペクトログラムを図2に示す(スペクトログラムを表示する部分までのmatlabコードを下に示す).スペクトログラムを見ると信号の開始時点及び終了時点において,周波数の広い範囲に大きい値が出ている(図3).これは信号が突然0→1に立ち上がる(終了時点では1→0になる)ことから,その時点においてあらゆる周波数が反応するのは当然である.しかしA[Hz]の正弦信号を入力したらA[Hz]のところにだけ値が出て来て欲しいわけで,そういうときは普通ならばスペクトログラムのウィンドウ幅を長くするなど分析側のパラメータを調整する必要がある.今回は理想的な信号が入ってくる場合の分析器の挙動を見たいので,分析側のパラメータを変えずに,実験用の信号を理想的なものに調整して綺麗な結果を出してみる.

f:id:Yuri-Processing-Club:20181003162409p:plain:w400
図 1: 信号の時間-振幅グラフ

f:id:Yuri-Processing-Club:20181003161738p:plain:w400
図 2: 信号のスペクトログラム

f:id:Yuri-Processing-Club:20181003163917p:plain:w400
図 3: 信号のスペクトログラム(説明付)



% 
% 信号を生成してスペクトログラムを表示する
% 
clear all;
close all;
clc

Fs = 32000;
SNRdB = -10;

% 信号を生成する
t1 = 1.5;
f1 = 5000;
N = Fs * t1;
t = ((1:N)-1)/Fs;
s1 = sin(2*pi*f1*t);
t2 = 1.5;
f2 = 10000;
N = Fs * t2;
t = ((1:N)-1)/Fs;
s2 = sin(2*pi*f2*t);
s0 = horzcat(s1, s2);
% 信号の前後にpad_time秒間の無音部分をくっつける
pad_time = 1;
s = padarray(s0, [0 round(pad_time * Fs)], 0, 'both');

% 時間グラフを描く
N = length(s);
t = ((1:N)-1)/Fs;
figure
plot(t, s);
title('信号の時間-振幅グラフ');
xlabel('時間(sec)');
ylabel('振幅');

% 信号のスペクトログラムを描画
Nw = 1024;
Nol = Nw*3/4;
NFFT = Nw;
figure;
spectrogram(s, Nw, Nol, NFFT, Fs, 'yaxis'); %ハミング窓
title('信号のスペクトログラム');
xlabel('時間(sec)');
colormap 'jet'
caxis_max = -20;
caxis_min = -120;
caxis([caxis_min, caxis_max]);



つまりこういうことをしたいわけである.

f:id:Yuri-Processing-Club:20181003165529p:plain:w500
図 4: 立ち上がり立ち下りが滑らかな信号の時間-振幅グラフ

信号の最初と終わりの部分だけ窓関数を掛けたような形である.ネタバレするとコサインテーパーウィンドウを掛けてやればいいわけだが,最初はハンウィンドウを前半後半に分けて信号の最初と終わりにそれぞれ掛けるというやり方を取っていた.

f:id:Yuri-Processing-Club:20181003171032p:plain:w400
図 5: ハンウィンドウ

f:id:Yuri-Processing-Club:20181003171715p:plain
図 6: ハンウィンドウを分割して信号の前後とかけあわせるやつの説明

%
% 上のコードの「信号を生成する」 の部分を下に入れ替える
%
% 信号の両端に少しだけ窓関数をかける
single_sided_window_length_sec = 0.2;% 片側窓関数がかけられる時間
single_sided_window_length = single_sided_window_length_sec * Fs;% 窓関数の片側の長さ
w = hann(single_sided_window_length * 2);% ハン窓を作成
figure
N = length(w);
t = ((1:N)-1)/Fs;
plot(t, w);
title('ハンウィンドウ');
xlabel('時間(sec)');

% 信号を生成する
t1 = 1.5;
f1 = 5000;
N = Fs * t1;
t = ((1:N)-1)/Fs;
s1 = sin(2*pi*f1*t);
% 窓関数をかける
s1(1:single_sided_window_length) = w(1:single_sided_window_length)' .* s1(1:single_sided_window_length);% 左側
s1(end-single_sided_window_length+1:end) = w(single_sided_window_length+1:end)' .* s1(end-single_sided_window_length+1:end);% 右側

t2 = 1.5;
f2 = 10000;
N = Fs * t2;
t = ((1:N)-1)/Fs;
s2 = sin(2*pi*f2*t);
% 窓関数をかける
s2(1:single_sided_window_length) = w(1:single_sided_window_length)' .* s2(1:single_sided_window_length);% 左側
s2(end-single_sided_window_length+1:end) = w(single_sided_window_length+1:end)' .* s2(end-single_sided_window_length+1:end);% 右側

s0 = horzcat(s1, s2);

% 信号の前後にpad_time秒間の無音部分をくっつける
pad_time = 1;
s = padarray(s0, [0 round(pad_time * Fs)], 0, 'both');

ハンウィンドウを分割して信号の前後にかけたものの時間-振幅グラフを図7,スペクトログラムを図8にそれぞれ示す.0→1に滑らかに変化しているので,スペクトログラムで欲しいところだけに信号が表示されている.

f:id:Yuri-Processing-Club:20181003174159p:plain:w400
図 7: 立ち上がり立ち下りを滑らかにした信号の時間-振幅グラフ

f:id:Yuri-Processing-Club:20181003174145p:plain:w400
図 8: 立ち上がり立ち下りを滑らかにした信号のスペクトログラム



こんなめんどくさいこと,やる?と思い上司に上手いやり方は無いか相談したところ,そういう場合は普通はコサインテーパーを使うとのこと.jp.mathworks.com
コサインテーパーウィンドウは開始部分及び終了部分がコサインの形をしていて,真ん中あたりは値が1のウィンドウである.信号と同じ長さのコサインテーパーウィンドウを用意して,信号にかけてやればOKである.
f:id:Yuri-Processing-Club:20181003180521p:plain:w400
図 9: 立ち上がり立ち下りを滑らかにした信号の時間-振幅グラフ(コサインテーパー使用)

%
% 上のコードの「信号を生成する」 の部分を下に入れ替える
%
single_sided_window_length_sec = 0.2;% 片側窓関数がかけられる時間

% 信号を生成する
t1 = 1.5;
f1 = 5000;
N = Fs * t1;
t = ((1:N)-1)/Fs;
s1 = sin(2*pi*f1*t);
% 窓関数をかける
window_signal_length_ratio = (single_sided_window_length_sec*2)/t1;% 窓関数がかけられる時間の比
w = tukeywin(N, window_signal_length_ratio);
figure
plot(t, w);
title('コサインテーパーウィンドウ');
xlabel('時間(sec)');
s1 = w' .* s1;

t2 = 1.5;
f2 = 10000;
N = Fs * t2;
t = ((1:N)-1)/Fs;
s2 = sin(2*pi*f2*t);
% 窓関数をかける
window_signal_length_ratio = (single_sided_window_length_sec*2)/t2;% 窓関数がかけられる時間の比
w = tukeywin(N, window_signal_length_ratio);
% 5kHzのコサインテーパーウィンドウと同じなのでグラフは省略
s2 = w' .* s2;

s0 = horzcat(s1, s2);

% 信号の前後にpad_time秒間の無音部分をくっつける
pad_time = 1;
s = padarray(s0, [0 round(pad_time * Fs)], 0, 'both');

f:id:Yuri-Processing-Club:20181003180010p:plain:w400
図 10: コサインテーパーウィンドウ

f:id:Yuri-Processing-Club:20181003175920p:plain:w400
図 11: 立ち上がり立ち下りを滑らかにした信号の時間-振幅グラフ(コサインテーパー使用)

f:id:Yuri-Processing-Club:20181003180853p:plain:w400
図 12: 立ち上がり立ち下りを滑らかにした信号のスペクトログラム(コサインテーパー使用)