窓関数

参考文献及びWEB

参考文献

(1)「ディジタル信号処理技術」玉井徳みち、長島厚、藤田泰弘、若井修造著日経BP社
(2)「ディジタル信号処理の基礎」三上直樹著 CQ出版
(3)「C言語によるディジタル信号処理入門」三上直樹著 CQ出版
(4)「アナログ&ディジタルフィルタ入門」小野浩司著日刊工業
(5)「フーリエの冒険」ヒッポファミリークラブ

あわせて読む

窓関数とは？

上図は、サンプリング周期 0.01[s]で5[Hz]の正弦波を
左は5周期ぴったり、右は5.5周期出力した図である。

フーリエ変換した結果、
左は5[Hz]で出力しているが、右は5Hz付近を中心にスペクトルが漏れている。

フーリエ変換では、信号の切り出し区間をうまく選べば、スペクトルの漏れを
なくすことができるが、一般に信号は多くの周波数成分を含むため切り出し区間をうまく選ぶことは困難である。

Matlab コード

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab sin spectram sample program
%                 programming by embedded.samurai
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

%サンプリング周期 0.01[s]
dt=0.01;

%サンプリング周波数100Hz 50Hz以下まで表すことができる
df=1/dt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1Hz正弦波を表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,1);

f1=1;
f2=2;
f3=5; %5[Hz]

T = 1/f3;
T_num = T/dt

% N=T_num(1周期の個数) x 個数分
N=T_num * 5
n=[0:1:N-1];

y1=100*sin(2*pi*f1*n*dt);
y2=100*sin(2*pi*f2*n*dt);
y3=80*sin(2*pi*f3*n*dt);

%y=y1+y2+y3;
y=y3;

t=n*dt;

% 1Hzの正弦波を表示する
plot(t,y,'r-','linewidth',2);

axis([0 Inf -Inf Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Time[s]','Fontsize',15);
ylabel('Amplitude','Fontsize',15)

y_w = y;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 周波数スペクトルを表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,3);

%基本周波数の導出 基本周波数=サンプリング周波数/データ数
basicfre = df/N

%FFTの結果は複素数である。複素数はテキストデータで保存できない
fftdata = fft(y_w(1:N))';

%離散周波数を作る
fre = [0:N-1]'*df/N;

%FFTの絶対値を求める
power=abs(fftdata).^2;
maxpower = max(power);
powerDB=20*log10(power/maxpower);
  
stem(fre,power,'r-*','linewidth',2);

axis([0 10  0 Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Frequency[Hz]','Fontsize',15);
ylabel('Normalize','Fontsize',15);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1Hz正弦波を表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,2);

f3=5; %5[Hz]

T = 1/f3;
T_num = T/dt

% N=T_num(1周期) x 個数分
N= T_num * 5.5;
n=[0:1:N-1];

y1=100*sin(2*pi*f1*n*dt);
y2=100*sin(2*pi*f2*n*dt);
y3=80*sin(2*pi*f3*n*dt);

%y=y1+y2+y3;
y=y3;

t=n*dt;

% 1Hzの正弦波を表示する
plot(t,y,'b-','linewidth',2);
axis([0 Inf -Inf Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Time[s]','Fontsize',15);
ylabel('Amplitude','Fontsize',15)

y_w = y;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 周波数スペクトルを表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,4);

%基本周波数の導出 基本周波数=サンプリング周波数/データ数
basicfre = df/N

%FFTの結果は複素数である。複素数はテキストデータで保存できない
fftdata = fft(y_w(1:N))';

%離散周波数を作る
fre = [0:N-1]'*df/N;

%FFTの絶対値を求める
power=abs(fftdata).^2;
maxpower = max(power);
powerDB=20*log10(power/maxpower);

stem(fre,power,'b-*','linewidth',2);

axis([0 10  0 Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Frequency[Hz]','Fontsize',15);
ylabel('Normalize','Fontsize',15);

print -djpeg spectram_sin.jpg

スペクトルの漏れを軽減するためには？

スペクトルの漏れを軽減するために窓関数と呼ばれるものが考案された。

このような窓関数とFFTする元信号を掛け合わせることにより、
スペクトルの漏れを軽減する。

今、元信号を s(n) とし、点で切り出したとき
窓関数の範囲は、 $0 \leq n \leq N-1$
となる。

窓関数は上から、

(a) 矩形窓
何も窓関数を使わないときは、矩形窓をかけたものとみなすことができる。
(元信号データの両端が突然0になることと同じだからである)

w(n) = 1

(b) ハニング窓
$w(n) = 0.5-0.5\cos\frac{2\pi n}{N-1}$

(d) ブラックマン窓
$w(n) = 0.42 - 0.5\cos\frac{2\pi n}{N-1} + 0.08\cos\frac{4\pi n}{N-1}$

実際に窓関数を適用する

さきほどの5.5周期の5Hz正弦波を

(1)矩形窓をとおした場合（つまり何も通さない場合)を赤
w(n) = 1

(2)ハニング窓をとおした場合を青
$w(n) = 0.5-0.5\cos\frac{2\pi n}{N-1}$

とします。
実際に窓関数をかけるには、

$s_{w}(n) = s(n) \times w(n)$

とします(Matlabソースも参考にしてください)。

このように窓関数をかけると周波数の漏れが抑制されているのが分かります。

窓関数の特徴

	優<--------------------->劣
周波数分解能	ハミングハニングブラックマン
スペクトルの漏れの抑制	ブラックマンハニングハミング

用例: $s(t) = A \sin(2\pi f_{A} t )+ B \sin(2\pi f_{B} t )$

(1) f_A と f_B が近い(a) AとBが同じくらい -----> 分解能を優先(b) AとBが大きく異なる -----> f_A と f_B を区別することは難しい

(2) f_A と f_B が遠い(a) AとBが同じくらい -----> どの窓関数でも良い(b) AとBが大きく異なる -----> 漏れの抑制を優先

Matlabソース

窓関数一覧

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab Window function
%                 programming by embedded.samurai
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

%サンプリング周期 0.01[s]
dt=0.01;

%サンプリング周波数100Hz 50Hz以下まで表すことができる
df=1/dt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 窓関数を表示する 矩形窓
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(4,1,1);

f3=5; %5[Hz]

T = 1/f3;
T_num = T/dt

% T_num(1周期の個数) x 個数分
N=T_num * 5
n=[0:1:N-1];

data=[];
for x=n
  tmp=1;
  data=[data,tmp];
end

%擬似矩形窓(くけいまど)を作る
plot(n,data,'r-','linewidth',8);
hold on;

for x=[0:0.01:1]
 plot(0,x,'rs-','linewidth',3,...
             'MarkerEdgeColor','r',...
             'MarkerFaceColor','r',...
             'MarkerSize',2);
 hold on
 plot(N-1,x,'rs-','linewidth',3,...
             'MarkerEdgeColor','r',...
             'MarkerFaceColor','r',...
             'MarkerSize',2);
 hold on
end

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

axis([0 N 0 1]);

xlabel('Number','Fontsize',15);
ylabel('Amplitude','Fontsize',15)
title('矩形窓');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 窓関数を表示する ハニング窓
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(4,1,2);

hanning  =0.5  -  0.5*cos(2*pi*n/(N-1));
hamming  =0.54 - 0.46*cos(2*pi*n/(N-1));
blackman =0.42 - 0.5*cos(2*pi*n/(N-1)) + 0.08*cos(4*pi*n/(N-1));

window_func = hanning;

plot(n,window_func,'r-','linewidth',3);
              
h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

axis([0 N 0 1]);

xlabel('Number','Fontsize',15);
ylabel('Amplitude','Fontsize',15)
title('ハニング窓');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 窓関数を表示する ハミング窓
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(4,1,3);

hanning  =0.5  -  0.5*cos(2*pi*n/(N-1));
hamming  =0.54 - 0.46*cos(2*pi*n/(N-1));
blackman =0.42 - 0.5*cos(2*pi*n/(N-1)) + 0.08*cos(4*pi*n/(N-1));

window_func = hamming;

plot(n,window_func,'r-','linewidth',3);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

axis([0 N 0 1]);

xlabel('Number','Fontsize',15);
ylabel('Amplitude','Fontsize',15);
title('ハミング窓');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 窓関数を表示する ブラックマン窓
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(4,1,4);

hanning  =0.5  -  0.5*cos(2*pi*n/(N-1));
hamming  =0.54 - 0.46*cos(2*pi*n/(N-1));
blackman =0.42 - 0.5*cos(2*pi*n/(N-1)) + 0.08*cos(4*pi*n/(N-1));

window_func = blackman;

plot(n,window_func,'r-','linewidth',3);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

axis([0 N 0 1]);

xlabel('Number','Fontsize',15);
ylabel('Amplitude','Fontsize',15)
title('ブラックマン窓');

print -djpeg window_function_gaiyou.jpg

5Hzの正弦波を窓関数をかけてフーリエ変換する wfunction_sin.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matlab FFTxWindow Function  sample program
%                 programming by embedded.samurai
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

%サンプリング周期 0.01[s]
dt=0.01;

%サンプリング周波数100Hz 50Hz以下まで表すことができる
df=1/dt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1Hz正弦波を表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(4,2,1);

f3=5; %5[Hz]

T = 1/f3;
T_num = T/dt

% T_num(1周期の個数) x 個数分
N=T_num * 5.5
n=[0:1:N-1];

y3=100*sin(2*pi*f3*n*dt);
y=y3;

t=n*dt;

% 1Hzの正弦波を表示する
plot(t,y,'r-','linewidth',3);

axis([0 Inf -Inf Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Time[s]','Fontsize',15);
ylabel('Amplitude','Fontsize',15)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 窓関数を表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(4,2,3);
n=[0:1:N-1];

data=[];
for x=n
  tmp=1;
  data=[data,tmp];
end

%擬似矩形窓(くけいまど)を作る
plot(n,data,'r-','linewidth',8);
hold on;

for x=[0:0.01:1]
 plot(0,x,'rs-','linewidth',3,...
              'MarkerEdgeColor','r',...
              'MarkerFaceColor','r',...
              'MarkerSize',2);
 hold on
 plot(N-1,x,'rs-','linewidth',3,...
             'MarkerEdgeColor','r',...
             'MarkerFaceColor','r',...
             'MarkerSize',2);
 hold on
end

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

axis([0 N 0 1]);

xlabel('Number','Fontsize',15);
ylabel('Amplitude','Fontsize',15)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 窓関数x正弦波を表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(4,2,5);

y_w = y; 

% 1Hzの正弦波を表示する
plot(t,y_w,'r-','linewidth',3);

axis([0 Inf -Inf Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Time[s]','Fontsize',15);
ylabel('Amplitude','Fontsize',15)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 周波数スペクトルを表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(4,2,7);

%基本周波数の導出 基本周波数=サンプリング周波数/データ数
basicfre = df/N

%FFTの結果は複素数である。複素数はテキストデータで保存できない
fftdata = fft(y_w(1:N))';

%離散周波数を作る
fre = [0:N-1]'*df/N;

%FFTの絶対値を求める
power=abs(fftdata).^2;
maxpower = max(power);
powerDB=20*log10(power/maxpower);

stem(fre,power,'r-*','linewidth',3);
 axis([0 10  0 Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Frequency[Hz]','Fontsize',15);
ylabel('Normalize','Fontsize',15);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1Hz正弦波を表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(4,2,2);

f3=5; %5[Hz]

T = 1/f3;
T_num = T/dt

% T_num(1周期の個数 20) x 個数分
n=[0:1:N-1];

y3=80*sin(2*pi*f3*n*dt);
y=y3;

t=n*dt;

% 1Hzの正弦波を表示する
plot(t,y,'b-','linewidth',3); 
axis([0 Inf -Inf Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Time[s]','Fontsize',15);
ylabel('Amplitude','Fontsize',15)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 窓関数を表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(4,2,4);

hanning  =0.5  -  0.5*cos(2*pi*n/(N-1));
hamming  =0.54 - 0.46*cos(2*pi*n/(N-1));
blackman =0.42 - 0.5*cos(2*pi*n/(N-1)) + 0.08*cos(4*pi*n/(N-1));

window_func = hanning;

plot(n,window_func,'b-','linewidth',3);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

axis([0 N -Inf Inf]);

xlabel('Number','Fontsize',15);
ylabel('Amplitude','Fontsize',15)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 窓関数x正弦波を表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(4,2,6); 

y_w = y.*window_func;

% 1Hzの正弦波を表示する
plot(t,y_w,'b-','linewidth',3);

axis([0 Inf -Inf Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Time[s]','Fontsize',15);
ylabel('Amplitude','Fontsize',15)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 周波数スペクトルを表示する
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(4,2,8);

%基本周波数の導出 基本周波数=サンプリング周波数/データ数
basicfre = df/N

%FFTの結果は複素数である。複素数はテキストデータで保存できない
fftdata = fft(y_w(1:N))';

%離散周波数を作る
fre = [0:N-1]'*df/N;

%FFTの絶対値を求める
power=abs(fftdata).^2;
maxpower = max(power);
powerDB=20*log10(power/maxpower);

stem(fre,power,'b-*','linewidth',3);
axis([0 10  0 Inf]);

h=gca;
set(h,'LineWidth',2,...,
'FontSize',15);

xlabel('Frequency[Hz]','Fontsize',15);
ylabel('Normalize','Fontsize',15);
  
print -djpeg window_function_sin.jpg