怎样在MATLAB中实现这几种常见的窗函数

uperrua

怎样在MATLAB中实现这几种常见的窗函数
" f2 h( j1 v4 w4 n
3 o: J9 e' E$ z4 i6 Y
+ |. P  M. x0 c5 P) S数字信号处理中通常是取其有限的时间片段进行分析，而不是对无限长的信号进行测量和运算。具体做法是从信号中截取一个时间片段，然后对信号进行傅里叶变换、相关分析等数学处理。信号的截断产生了能量泄漏，而用FFT算法计算频谱又产生了栅栏效应，从原理上讲这两种误差都是不能消除的。在FFT分析中为了减少或消除频谱能量泄漏及栅栏效应，可采用不同的截取函数对信号进行截短，截短函数称为窗函数，简称为窗。
( a1 d8 f) ?( U3 a( u; A泄漏与窗函数频谱的两侧旁瓣有关，对于窗函数的选用总的原则是，要从保持最大信息和消除旁瓣的综合效果出发来考虑问题，尽可能使窗函数频谱中的主瓣宽度应尽量窄，以获得较陡的过渡带；旁瓣衰减应尽量大，以提高阻带的衰减，但通常都不能同时满足这两个要求。频谱中的如果两侧瓣的高度趋于零，而使能量相对集中在主瓣，就可以较为接近于真实的频谱。不同的窗函数对信号频谱的影响是不一样的，这主要是因为不同的窗函数，产生泄漏的大小不一样，频率分辨能力也不一样。信号的加窗处理，重要的问题是在于根据信号的性质和研究目的来选用窗函数。图1是几种常用的窗函数的时域和频域波形，其中矩形窗主瓣窄，旁瓣大，频率识别精度最高，幅值识别精度最低，如果仅要求精确读出主瓣频率，而不考虑幅值精度，则可选用矩形窗，例如测量物体的自振频率等；布莱克曼窗主瓣宽，旁瓣小，频率识别精度最低，但幅值识别精度最高；如果分析窄带信号，且有较强的干扰噪声，则应选用旁瓣幅度小的窗函数，如汉宁窗、三角窗等；对于随时间按指数衰减的函数，可采用指数窗来提高信噪比。表1 是几种常用的窗函数的比较。
$ A! e8 y1 B- W6 Q; J4 d- G
如果被测信号是随机或者未知的，或者是一般使用者对窗函数不大了解，要求也不是特别高时，可以选择汉宁窗，因为它的泄漏、波动都较小，并且选择性也较高。但在用于校准时选用平顶窗较好，因为它的通带波动非常小，幅度误差也较小。
1 J! L. c/ _, M& R+ }. ?% I, j
7 f7 Y& B9 i" f( q  m

表1 几种常用的窗函数的比较

+ [3 O$ {' L5 c3 w

名称	特点	应用
矩形窗 Rectangle	矩形窗使用最多，习惯上不加窗就是使信号通过了矩形窗。这种窗的优点是主瓣比较集中，缺点是旁瓣较高，并有负旁瓣，导致变换中带进了高频干扰和泄漏，甚至出现负谱现象。频率识别精度最高，幅值识别精度最低，所以矩形窗不是一个理想的窗。	如果仅要求精确读出主瓣频率，而不考虑幅值精度，则可选用矩形窗，例如测量物体的自振频率等，也可以用在阶次分析中。
汉宁窗 Hanning	又称升余弦窗。主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗．但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨力下降。它与矩形窗相比，泄漏、波动都减小了,并且选择性也提高。	是很有用的窗函数。如果测试信号有多个频率分量，频谱表现的十分复杂，且测试的目的更多关注频率点而非能量的大小，需要选择汉宁窗。如果被测信号是随机或者未知的，选择汉宁窗。
海明窗 （汉明窗） Hamming	与汉宁窗都是余弦窗，又称改进的升余弦窗，只是加权系数不同，使旁瓣达到更小。但其旁瓣衰减速度比汉宁窗衰减速度慢。	与汉明窗类似，也是很有用的窗函数。
平顶窗 Flap Top	平顶窗在频域时的表现就象它的名称一样有非常小的通带波动。	由于在幅度上有较小的误差，所以这个窗可以用在校准上。
凯塞窗 Kaiser	定义了一组可调的由零阶贝塞尔Bessel 函数构成的窗函数，通过调整参数β可以在主瓣宽度和旁瓣衰减之间自由选择它们的比重。对于某一长度的Kaiser 窗，给定β，则旁瓣高度也就固定了。
布莱克曼窗 Blackman	二阶升余弦窗，主瓣宽，旁瓣比较低，但等效噪声带宽比汉宁窗要大一点，波动却小一点。频率识别精度最低，但幅值识别精度最高，有更好的选择性。	常用来检测两个频率相近幅度不同的信号。
高斯窗 Gaussian	是一种指数窗。主瓣较宽，故而频率分辨力低；无负的旁瓣，第一旁瓣衰减达一55dB。常被用来截短一些非周期信号，如指数衰减信号等。	对于随时间按指数衰减的函数，可采用指数窗来提高信噪比。
三角窗 （费杰窗） Fejer	是幂窗的一次方形式。与矩形窗比较，主瓣宽约等于矩形窗的两倍，但旁瓣小，而且无负旁瓣。	如果分析窄带信号，且有较强的干扰噪声，则应选用旁瓣幅度小的窗函数，如汉宁窗、三角窗等；
切比雪夫窗（Chebyshev）	在给定旁瓣高度下，Chebyshev窗的主瓣宽度最小，具有等波动性，也就是说，其所有的旁瓣都具有相等的高度。

2 G  v6 ~3 F2 V
: Z1 L' F  {7 w8 N' l7 x1 `. H
& \" x/ G8 y1 Z6 e下面是几种窗函数归一化DTFT幅度的MATLAB程序:
6 l% B% ], ~0 q  f; c; O
3 N9 _$ T5 O/ s7 Z附上DTFT函数（dtft.m）：
1 ^( i& }% B: V# _5 q. F0 q( E/ w0 s
function [ X ] = dtft( x,n,w )

%   Computes Discrete-time Fourier Transform
4 z8 J  m- k" q- P$ U
9 X# k" j# a  w% o% a%   [X] = dtft(x,n,w)
0 S, {5 q8 t5 \+ b4 V
%    X = DTFT values computed at w.frequencies
& i; T5 J# t6 x- e! j" ]
%    x = finite duration sequence over n

%    n = sample position vector
+ B: x  e2 e& ]) c
%    w = frequency location vector
4 u  R; a! U$ Z3 V! ^
. Y; v7 d* {3 pX = x*exp(-j*n'*w);
+ |( X" B- s! O, ^, k/ @8 _
%
5 v# I  r; x$ t9 h
( g: W2 h: Z) d: T1 f3 Nend


) a! Z9 V0 W4 v$ o8 ~1 K
矩形窗:

%DTFT of a Rectangular Window, M=10,25,50,101

( A. j5 F; n4 N2 t/ n, R! v: l8 [clc; close all;
$ n2 a9 R% K% Q6 c2 n% @* R
) m9 I% w. h7 ?) S# ~& z: t; y2 NHf_1=figure; set(Hf_1,'NumberTitle','off','Name','P0304a');
( s- o3 g; R+ L/ u9 R
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));

, _- j1 J! n/ U5 v9 K4 vX=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]);

; [3 x6 [) I3 h8 N4 _ylabel('|X|'); title(['M=10']);
! Z5 t- [" P% a# B6 V! L; }
set(gca,'XTick',wtick,'YTick',magtick);
$ J& @: p$ q, i
4 \9 `! r# q0 N, V' S) H2 @/ t% M=25
0 f- a; K  h- V/ Q6 m
M=25; n=0:M; x=ones(1,length(n));

X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
! i: H4 K. }8 z' p) p
" W( }5 Q; F" \. b" F8 Psubplot(2,2,2); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
" t/ f' }! U$ @- [
title(['M=25']); set(gca,'XTick',wtick,'YTick',magtick);
  U5 k4 }$ ]1 b
1 |( u5 J  l  X5 X! z% M=50

" p; [$ \8 g5 c, [/ kM=50; n=0:M; x=ones(1,length(n));
3 I  O- ^6 k  Z. S9 p6 c
1 F+ y1 Z, w  a4 F( e5 ?X=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
: {6 O$ C7 [; x" b# C9 O9 U; @
' m) v3 u) ^, zsubplot(2,2,3); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);

" u; F, ]% D) Exlabel('\omega/\pi'); ylabel('|X|');

title('M=50'); set(gca,'XTick',wtick,'YTick',magtick);

& y5 z, V& a! }: N2 N$ v& Q% M=101

0 K: v, t% {% M, A9 D) U& N1 G- xM=101; n=0:M; x=ones(1,length(n));
$ C- l% V' z0 K
) J" O- ~+ W$ n5 P7 AX=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
7 @6 d! M6 d' D. w
subplot(2,2,4); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
1 p( u* S, }. D. g
xlabel('\omega/\pi'); ylabel('|X|');
. p( s  n% A0 v# g5 [; N- ]$ y
  a( e4 g7 y+ R0 ftitle(['M=101']);


4 r+ B3 ?, V; w三角窗:

% Triangular Window:

1 `3 {4 v$ |: u* k$ v  V; {% DTFT of a Triangular Window,M = 10,25,50,101

clc; close all;

Hf_1=figure; set(Hf_1,'NumberTitle','off','Name','P0304b');

2 w& o/ U3 V6 ow=linspace(-pi,pi,501); wtick=[-1:0.5:1]; magtick=[0:0.5:1.1];
0 h5 A% }8 E, M7 ]* I
% M = 10
$ q, A4 H9 ]- X: w
0 o" x7 k; v2 j6 ^) ]M=10; n=0:M; x=(1-(abs(M-1-(2*n))/(M+1)));
  R6 g( U* }# p" n8 d, D" ~
3 n+ L  [/ \  B1 nx=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);

+ g# i$ T9 n' Y9 Asubplot(2,2,1); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
4 I! L1 N& E6 c& f3 \* c/ W! ]9 e
ylabel('|X|'); title(['M = 10']);

" P  O, H0 I; S+ h4 k& g8 b7 |. Pset(gca,'XTick',wtick,'YTick',magtick);

% M = 25
' F& |" E. @; M
/ ^! |/ ]1 g- SM=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]);

3 Y' h$ I! ]6 `6 H; K6 Qtitle(['M = 25']); set(gca,'XTick',wtick,'YTick',magtick);

% M = 50
+ b9 ~1 F+ [% I
0 D4 R9 [+ z; z9 ~2 LM=50; n=0:M; x=(1-(abs(M-1-(2*n))/(M+1)));
  m' W( s1 V  h
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]);

( n6 L8 i# u* Y  ^xlabel('\omega/\pi'); ylabel('|X|'); title(['M = 50']);
8 n1 h" F9 E, b2 t: v& F
! @* s/ _4 N* D$ y: L4 t# s; wset(gca,'XTick',wtick,'YTick',magtick);
9 r$ Y$ t6 M/ i0 J  F, J7 Y
% M = 100
6 G) W- A* c: I$ n( }
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);
+ D0 y( y% T( S% J8 v- z
+ `/ U0 `$ b' g+ ]2 E4 Tsubplot(2,2,4); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);

- R4 |$ O# d% ?: S5 {xlabel('\omega/\pi'); title(['M = 101']);

set(gca,'XTick',wtick,'YTick',magtick);
0 x" i) \# G# H1 _1 |( b7 p- {

& I" N: @& J' ~' {: y  X
! ]+ O* f3 ~5 H+ d$ y! {海宁窗:

% Hann Window

3 a: f. h  q, P+ T% DTFT of a Hann Window, M = 10,25,50,101
) w3 V( H% @8 m! R$ u1 D
2 d+ G% ^+ I  [8 B% e; P& i3 yclc;close all;

; O, v  u$ ?2 s* b; Y  ^8 b# jHf_1 = figure; set(Hf_1,'NumberTitle','off','Name','P0304c');

6 j" j; h: t& @w=linspace(-pi,pi,501); wtick=[-1:0.5:1]; magtick=[0:0.5:1.1];
+ \7 Q, A& x# B; f' e4 e( x
% M = 10

3 `' V) H( A# J- l+ M4 T% ?& kM=10;n=0:M; x=0.5*(1-cos((2*pi*n)/(M-1)));
9 }7 J) `, w- A$ i* f
1 V7 T5 v( m, I/ yX=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]);
& ^. v) V# K, z* c7 ]1 h: w+ J
4 O6 ?6 _; u2 Q* W' }) Tylabel('|X|'); title(['M = 10']);
( s! p9 m: b1 f$ g
2 z! R4 N0 i6 k5 Cset(gca,'XTick',wtick,'YTick',magtick);

7 b9 G4 e3 }+ @3 R/ g% M = 25

M=25;n=0:M; x=0.5*(1-cos((2*pi*n)/(M-1)));
) {2 N. d' V/ h* ]8 b9 ]' o2 l
* J/ q, a1 K! z, n% W9 i7 jX=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);
- J/ ^* y! w, Q* J0 }
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);

+ g0 |3 o" z. v  w2 z% M = 50

7 ?, @; R/ R. o# `/ [- UM=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]);
& W% p9 X6 @/ I+ r, p
xlabel('\omega/\pi'); ylabel('|X|'); title(['M = 50']);

% O" b* a* E8 s+ d/ Dset(gca,'XTick',wtick,'YTick',magtick);

% M = 101

& m% x! f+ z; L3 M+ fM=101;n=0:M; x=0.5*(1-cos((2*pi*n)/(M-1)));
* Y, n0 ?* K: a  D2 r1 V
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]);

& T  h; G- G  ^, qxlabel('\omega/\pi'); title(['M = 101']);
* v% A0 \9 l0 _( W# v0 m# A
- T) H4 S9 v+ _& G5 iset(gca,'XTick',wtick,'YTick',magtick);

3 Z/ x, S- r5 }+ G
哈明窗:
( V, _( x4 N8 ]" j
0 W. b: R! `% A# E% Hamming Window:
8 _: f3 d, c, P8 U8 x5 o3 p& r8 P
9 O0 J) T' C* F6 qclc; close all;
; A/ w! b: ]1 x6 V
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

  h: d& }9 p; m  JM=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);
  e3 r: [/ O, D* }& u7 J
subplot(2,2,1); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
& ?5 J5 T2 e$ M) J' @
8 U# o' |# s( ~* u: X, I ylabel('|X|'); title(['M = 10']);

8 k8 _- J; s$ A9 y$ Sset(gca,'XTick',wtick,'YTick',magtick);

/ C: t' E  R9 b) H: r% M = 25
1 h  m/ I6 c1 h1 ^7 U3 c& e. b
$ g/ z/ u+ y! |M=25; n=0:M; x=(0.54-0.46*cos((2*pi*n)/(M-1)));
9 P( `! ]% ~. u8 g) I. K
( a; U" `$ W3 EX=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);

6 Y# r( A" Z' J8 Z6 a+ F2 D* ^  Hsubplot(2,2,2); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);
9 f9 |: L% ~5 \: `- x( T2 n
' O6 _" \, B2 j  A, [title(['M = 25']);

set(gca,'XTick',wtick,'YTick',magtick);

% M = 50
, J! g7 j, T, r# S/ C& e
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]);

; m9 ?( c; E0 q# O' b1 O- Oxlabel('\omega/\pi'); ylabel('|X|'); title(['M=50']);
7 A# D, e; d! z4 e6 g
set(gca,'XTick',wtick,'YTick',magtick);
! {8 K) {8 k1 \. Z9 p* P8 |. U
% M = 101
2 z" N6 t" F  W& H* f+ t
M=101; n=0:M; x=(0.54-0.46*cos((2*pi*n)/(M-1)));
3 f' ?% n3 G3 x9 J9 p( A
6 Q3 _$ P/ z% ?; m! LX=dtft(x,n,w); magX=abs(X); magX=magX/max(magX);

$ t( b9 I& V  N0 a0 K1 ?subplot(2,2,4); plot(w/pi,magX,'LineWidth',1.5); axis([-1 1 0 1.1]);

4 M# `: C3 x; Y& nxlabel('\omega/\pi');title(['M=101']);
+ Z: {/ E5 k! Q& c7 D  z* e
set(gca,'XTick',wtick,'YTick',magtick);

2 x0 X7 t/ }  Y4 e: W+ G' ~& `

: a6 @4 K% U/ p! h5 E/ q7 G* Z

relchhiclty

谢谢分享