幾種常見窗函式及其MATLAB實現

方克明發表於2017-03-30
函式Matlab
對於窗函式一些教材和論壇都有好的定義和解釋說明。鑑於前人的基礎上，我也作了一些學習和研究。本文針對於不同的窗函式進行簡單的對比，並有MATLAB程式應用。

注：原文來源：http://blog.csdn.net/demo_deng/article/details/17351647       感謝博主

數字訊號處理中通常是取其有限的時間片段進行分析，而不是對無限長的訊號進行測量和運算。具體做法是從訊號中擷取一個時間片段，然後對訊號進行傅立葉變換、相關分析等數學處理。訊號的截斷產生了能量洩漏，而用FFT演算法計算頻譜又產生了柵欄效應，從原理上講這兩種誤差都是不能消除的。在FFT分析中為了減少或消除頻譜能量洩漏及柵欄效應，可採用不同的擷取函式對訊號進行截短，截短函式稱為窗函式，簡稱為窗。


洩漏與窗函式頻譜的兩側旁瓣有關，對於窗函式的選用總的原則是，要從保持最大資訊和消除旁瓣的綜合效果出發來考慮問題，儘可能使窗函式頻譜中的主瓣寬度應儘量窄，以獲得較陡的過渡帶；旁瓣衰減應儘量大，以提高阻帶的衰減，但通常都不能同時滿足這兩個要求。頻譜中的如果兩側瓣的高度趨於零，而使能量相對集中在主瓣，就可以較為接近於真實的頻譜。不同的窗函式對訊號頻譜的影響是不一樣的，這主要是因為不同的窗函式，產生洩漏的大小不一樣，頻率分辨能力也不一樣。訊號的加窗處理，重要的問題是在於根據訊號的性質和研究目的來選用窗函式。圖1是幾種常用的窗函式的時域和頻域波形，其中矩形窗主瓣窄，旁瓣大，頻率識別精度最高，幅值識別精度最低，如果僅要求精確讀出主瓣頻率，而不考慮幅值精度，則可選用矩形窗，例如測量物體的自振頻率等；布萊克曼窗主瓣寬，旁瓣小，頻率識別精度最低，但幅值識別精度最高；如果分析窄帶訊號，且有較強的干擾噪聲，則應選用旁瓣幅度小的窗函式，如漢寧窗、三角窗等；對於隨時間按指數衰減的函式，可採用指數窗來提高訊雜比。表1 是幾種常用的窗函式的比較。


如果被測訊號是隨機或者未知的，或者是一般使用者對窗函式不大瞭解，要求也不是特別高時，可以選擇漢寧窗，因為它的洩漏、波動都較小，並且選擇性也較高。但在用於校準時選用平頂窗較好，因為它的通帶波動非常小，幅度誤差也較小。

                                              

表1 幾種常用的窗函式的比較

名稱

特點

應用

矩形窗

Rectangle

矩形窗使用最多，習慣上不加窗就是使訊號通過了矩形窗。這種窗的優點是主瓣比較集中，缺點是旁瓣較高，並有負旁瓣，導致變換中帶進了高頻干擾和洩漏，甚至出現負譜現象。頻率識別精度最高，幅值識別精度最低，所以矩形窗不是一個理想的窗。

如果僅要求精確讀出主瓣頻率，而不考慮幅值精度，則可選用矩形窗，例如測量物體的自振頻率等，也可以用在階次分析中。

漢寧窗

Hanning

又稱升餘弦窗。主瓣加寬並降低，旁瓣則顯著減小，從減小洩漏觀點出發，漢寧窗優於矩形窗．但漢寧窗主瓣加寬，相當於分析頻寬加寬，頻率分辨力下降。它與矩形窗相比，洩漏、波動都減小了,並且選擇性也提高。

是很有用的窗函式。如果測試訊號有多個頻率分量，頻譜表現的十分複雜，且測試的目的更多關注頻率點而非能量的大小，需要選擇漢寧窗。如果被測訊號是隨機或者未知的，選擇漢寧窗。

海明窗

（漢明窗）

Hamming

與漢寧窗都是餘弦窗，又稱改進的升餘弦窗，只是加權係數不同，使旁瓣達到更小。但其旁瓣衰減速度比漢寧窗衰減速度慢。

與漢明窗類似，也是很有用的窗函式。

平頂窗

Flap Top

平頂窗在頻域時的表現就象它的名稱一樣有非常小的通帶波動。

由於在幅度上有較小的誤差，所以這個窗可以用在校準上。

凱塞窗

Kaiser

定義了一組可調的由零階貝塞爾Bessel函式構成的窗函式，通過調整引數β可以在主瓣寬度和旁瓣衰減之間自由選擇它們的比重。對於某一長度的Kaiser窗，給定β，則旁瓣高度也就固定了。

 

布萊克曼窗

Blackman

二階升餘弦窗，主瓣寬，旁瓣比較低，但等效噪聲頻寬比漢寧窗要大一點，波動卻小一點。頻率識別精度最低，但幅值識別精度最高，有更好的選擇性。

常用來檢測兩個頻率相近幅度不同的訊號。

高斯窗

Gaussian

是一種指數窗。主瓣較寬，故而頻率分辨力低；無負的旁瓣，第一旁瓣衰減達一55dB。常被用來截短一些非週期訊號，如指數衰減訊號等。

對於隨時間按指數衰減的函式，可採用指數窗來提高訊雜比。

三角窗

（費傑窗）

Fejer

是冪窗的一次方形式。與矩形窗比較，主瓣寬約等於矩形窗的兩倍，但旁瓣小，而且無負旁瓣。

如果分析窄帶訊號，且有較強的干擾噪聲，則應選用旁瓣幅度小的窗函式，如漢寧窗、三角窗等。

切比雪夫窗Chebyshev

在給定旁瓣高度下，Chebyshev窗的主瓣寬度最小，具有等波動性，也就是說，其所有的旁瓣都具有相等的高度。

 


下面是幾種窗函式歸一化DTFT幅度的MATLAB程式:

附上DTFT函式（dtft.m）：

function [ X ] = dtft( x,n,w )
%   Computes Discrete-time Fourier Transform
%   [X] = dtft(x,n,w)
%    X = DTFT values computed at w.frequencies
%    x = finite duration sequence over n
%    n = sample position vector
%    w = frequency location vector
X = x*exp(-j*n'*w);
%
end


矩形窗:

%DTFT of a Rectangular Window, M=10,25,50,101
clc; close all;
Hf_1=figure; set(Hf_1,'NumberTitle','off','Name','P0304a');
w=linspace(-pi,pi,501); wtick=[-1:0.5:1]; magtick=[0:0.5:1.1];
% M=10
M=10; n=0:M; x=ones(1,length(n));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,1); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
ylabel('|X|'); title(['M=10']);
set(gca,'XTick',wtick,'YTick',magtick);
% M=25
M=25; n=0:M; x=ones(1,length(n));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,2); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
title(['M=25']); set(gca,'XTick',wtick,'YTick',magtick);
% M=50
M=50; n=0:M; x=ones(1,length(n));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,3); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
xlabel('\omega/\pi'); ylabel('|X|'); 
title('M=50'); set(gca,'XTick',wtick,'YTick',magtick);
% M=101
M=101; n=0:M; x=ones(1,length(n));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,4); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
xlabel('\omega/\pi'); ylabel('|X|');
title(['M=101']);

 

三角窗:

% Triangular Window: 
% DTFT of a Triangular Window,M = 10,25,50,101
clc; close all;
Hf_1=figure; set(Hf_1,'NumberTitle','off','Name','P0304b');
w=linspace(-pi,pi,501); wtick=[-1:0.5:1]; magtick=[0:0.5:1.1];
% M = 10
M=10; n=0:M; x=(1-(abs(M-1-(2*n))/(M+1)));
x=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,1); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
ylabel('|X|'); title(['M = 10']);
set(gca,'XTick',wtick,'YTick',magtick);
% M = 25
M=25; n=0:M; x=(1-(abs(M-1-(2*n))/(M+1)));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,2); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
title(['M = 25']); set(gca,'XTick',wtick,'YTick',magtick);
% M = 50
M=50; n=0:M; x=(1-(abs(M-1-(2*n))/(M+1)));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,3); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
xlabel('\omega/\pi'); ylabel('|X|'); title(['M = 50']); 
set(gca,'XTick',wtick,'YTick',magtick);
% M = 100
M=101;n=0:M; x=(1-(abs(M-1-(2*n))/(M+1)));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,4); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
xlabel('\omega/\pi'); title(['M = 101']); 
set(gca,'XTick',wtick,'YTick',magtick);

 

漢寧窗:

% Hann Window
% DTFT of a Hann Window, M = 10,25,50,101
clc;close all;
Hf_1 = figure; set(Hf_1,'NumberTitle','off','Name','P0304c');
w=linspace(-pi,pi,501); wtick=[-1:0.5:1]; magtick=[0:0.5:1.1];
% M = 10
M=10;n=0:M; x=0.5*(1-cos((2*pi*n)/(M-1)));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,1); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
ylabel('|X|'); title(['M = 10']);
set(gca,'XTick',wtick,'YTick',magtick);
% M = 25
M=25;n=0:M; x=0.5*(1-cos((2*pi*n)/(M-1)));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,2); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
 title(['M = 25']);
set(gca,'XTick',wtick,'YTick',magtick);
% M = 50
M=50;n=0:M; x=0.5*(1-cos((2*pi*n)/(M-1)));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,3); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
xlabel('\omega/\pi'); ylabel('|X|'); title(['M = 50']);
set(gca,'XTick',wtick,'YTick',magtick);
% M = 101
M=101;n=0:M; x=0.5*(1-cos((2*pi*n)/(M-1)));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,4); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
xlabel('\omega/\pi'); title(['M = 101']);
set(gca,'XTick',wtick,'YTick',magtick);


 

漢明窗:

% Hamming Window:
clc; close all;
Hf_1=figure; set(Hf_1,'NumberTitle','off','Name','P0304d');
w=linspace(-pi,pi,501); wtick=[-1:0.5:1]; magtick=[0:0.5:1.1];
% M = 10
M=10; n=0:M; x=(0.54-0.46*cos((2*pi*n)/(M-1)));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,1); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
 ylabel('|X|'); title(['M = 10']);
set(gca,'XTick',wtick,'YTick',magtick);
% M = 25
M=25; n=0:M; x=(0.54-0.46*cos((2*pi*n)/(M-1)));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,2); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
title(['M = 25']);
set(gca,'XTick',wtick,'YTick',magtick);
% M = 50
M=50; n=0:M; x=(0.54-0.46*cos((2*pi*n)/(M-1)));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,3); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
xlabel('\omega/\pi'); ylabel('|X|'); title(['M=50']);
set(gca,'XTick',wtick,'YTick',magtick);
% M = 101
M=101; n=0:M; x=(0.54-0.46*cos((2*pi*n)/(M-1)));
X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
subplot(2,2,4); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
xlabel('\omega/\pi');title(['M=101']);
set(gca,'XTick',wtick,'YTick',magtick);


相關文章