数字信号处理(MATLAB)[采样点数 频谱泄露 FIR滤波器设计]

文章目录

题目

如图1所示，由 x ( t ) x(t) x(t)信号 x 1 ( t ) x_1(t) x1​(t)、信号 x 2 ( t ) x_2(t) x2​(t)和信号 x 3 ( t ) x_3(t) x3​(t)的和构成。对 x ( t ) x(t) x(t)进行均匀采样（采样周期为 T s = 0.25 m s T_s=0.25ms Ts​=0.25ms），获得其样本信号 x ( n T s ) x(nT_s) x(nTs​)。

其中， x i ( t ) = A i c o s ( 2 π f i t ) x_i(t)=A_icos(2\pi f_it) xi​(t)=Ai​cos(2πfi​t)，参数 A i A_i Ai​和 f i f_i fi​如下表所示

问题

（供参考）

傅里叶变换

（1）求 x ( t ) x(t) x(t)的傅立叶变换，并画出其幅频响应曲线。

解：由于 x ( t ) x(t) x(t)由信号 x 1 ( t ) x_1(t) x1​(t)、信号 x 2 ( t ) x_2(t) x2​(t)和信号 x 3 ( t ) x_3(t) x3​(t)的和构成。由此可将 x ( t ) x(t) x(t)表示为:

x ( t ) = ∑ i = 1 3 x i ( t ) = ∑ i = 1 3 A i cos ⁡ ( 2 π f i t ) x(t)= \sum\limits_{i = 1}^3 {{x_i}\left( t \right)} {\rm{ = }}\sum\limits_{i = 1}^3 {{A_i}\cos ( {2\pi {f_i}t} )} x(t)=i=1∑3​xi​(t)=i=1∑3​Ai​cos(2πfi​t)

由傅里叶变换公式 X ( j ω ) = ∫ − ∞ + ∞ x ( t ) e − j ω t d t X(j\omega ) = \int_{ - \infty }^{ + \infty } {x(t)} {e^{ - j\omega t}}dt X(jω)=∫−∞+∞​x(t)e−jωtdt

F [ x ( t ) ] = F [ ∑ i = 1 3 x i ( t ) ] = ∑ i = 1 3 A i π [ δ ( ω − 2 π f i ) + δ ( ω + 2 π f i ) ] = 0.2 π [ δ ( ω − 2 ∗ 990 π ) + δ ( ω + 2 ∗ 990 π ) ] + 2 π [ δ ( ω − 2 ∗ 980 π ) + δ ( ω + 2 ∗ 980 π ) ] + 0.4 π [ δ ( ω − 2 ∗ 1600 π ) + δ ( ω + 2 ∗ 1600 π ) ] {\cal F}\left[ {x\left( t \right)} \right] = {\cal F}\left[ {\sum\limits_{i = 1}^3 {{x_i}\left( t \right)} } \right]\ = \sum\limits_{i = 1}^3 {{A_i}\pi \left[ {\delta (\omega - 2\pi {f_i}) + \delta (\omega + 2\pi {f_i})} \right]} \\ = 0.2\pi \left[ {\delta (\omega - 2*990\pi ) + \delta (\omega + 2*990\pi )} \right] + \\ 2\pi \left[ {\delta (\omega - 2*980\pi ) + \delta (\omega + 2*980\pi )} \right] + \\ 0.4\pi \left[ {\delta (\omega - 2*1600\pi ) + \delta (\omega + 2*1600\pi )} \right] F[x(t)]=F[i=1∑3​xi​(t)]=i=1∑3​Ai​π[δ(ω−2πfi​)+δ(ω+2πfi​)]=0.2π[δ(ω−2∗990π)+δ(ω+2∗990π)]+2π[δ(ω−2∗980π)+δ(ω+2∗980π)]+0.4π[δ(ω−2∗1600π)+δ(ω+2∗1600π)]

%傅里叶变换clc;clear;syms t;x1=0.2*cos(2*pi*1280*t);x2=2*cos(2*pi*1290*t);x3=0.4*cos(2*pi*1800*t);x=x1+x2+x3; %x(t)信号由x1、x2、x3的和组成。F=fourier(x)

%幅频响应曲线close all;clear;clc;Ts=0.00025; %采样周期，以0.25ms的周期均匀采样Fs=1/Ts; %采样频率t=0:Ts:1-Ts; %采样间隔时长x1=0.2*cos(2*pi*990*t); %三分量信号x2=2*cos(2*pi*980*t);x3=0.4*cos(2*pi*1600*t);x=x1+x2+x3; %x(t)信号由x1、x2、x3的和组成。figure(1);subplot(2,1,1);plot(t,x);xlabel('时间/s');ylabel('幅值');title("original signal");%绘制原始信号图%频域N=(1-Ts)/Ts+1;%采样点数n=0:N-1;fft_x=n*Fs/N-Fs/2;x_fft=fft(x);shift_f=abs(fftshift(x_fft));y2=shift_f/N;subplot(2,1,2);plot(fft_x,y2);xlabel('频率/Hz');ylabel('幅值');title('FFT'); %0频率分量位于坐标中心grid on;

x 3 ( t ) x_3(t) x3​(t)信号：

频谱图：

最少采样点数

(2)如果希望能从所获取到的样本信号 x ( n T s ) x(nT_s) x(nTs​)中，通过DFT有效的区分出信号 x 1 ( t ) x_1(t) x1​(t)、 x 2 ( t ) x_2(t) x2​(t)和 x 3 ( t ) x_3(t) x3​(t),试求最少需要的采样点数N，画出此时样本信号 x ( n T s ) x(nT_s) x(nTs​)（2）的离散傅里叶变换（DFT）的主值区间的频谱图，并标出数字频率及其对应的模拟频率。

解：分析可知原时域模拟信号为

x ( t ) = ∑ i = 1 3 x i ( t ) = ∑ i = 1 3 A i cos ⁡ ( 2 π f i t ) = 0.2 cos ⁡ ( 2 π × 990 t ) + 2 cos ⁡ ( 2 π × 980 t ) + 0.4 cos ⁡ ( 2 π × 1600 t ) x\left( t \right) = \sum\limits_{i = 1}^3 {{x_i}\left( t \right)} {\rm{ = }}\sum\limits_{i = 1}^3 {{A_i}\cos \left( {2\pi {f_i}t} \right)} = 0.2\cos \left( {2\pi \times 990t} \right) + 2\cos \left( {2\pi \times 980t} \right) + 0.4\cos \left( {2\pi \times 1600t} \right) x(t)=i=1∑3​xi​(t)=i=1∑3​Ai​cos(2πfi​t)=0.2cos(2π×990t)+2cos(2π×980t)+0.4cos(2π×1600t)

即信号的3个频率分量为 f 1 = 990 H z f_1=990Hz f1​=990Hz, f 2 = 980 H z f_2=980Hz f2​=980Hz, f 3 = 1600 H z f_3=1600Hz f3​=1600Hz

由采样周期为 T s = 0.25 m s T_s=0.25ms Ts​=0.25ms，根据抽样间隔与抽样频率的关系可知 f s = 1 T s = 1 0.25 m s = 4000 H z {f_s} = \frac{1}{{{T_s}}} = \frac{1}{{0.25ms}} = 4000Hz fs​=Ts​1​=0.25ms1​=4000Hz

抽样后得到的离散序列为 x ( n ) = x a ( t ) ∣ t = n T s = 0.2 cos ⁡ ( 99 π 200 n ) + 2 cos ⁡ ( 49 π 100 n ) + 0.4 cos ⁡ ( 4 π 5 n ) x(n) = {x_a}(t){|_{t = n{T_s}}} = 0.2\cos (\frac{{99\pi }}{{200}}n) + 2\cos (\frac{{49\pi }}{{100}}n) + 0.4\cos (\frac{{4\pi }}{5}n) x(n)=xa​(t)∣t=nTs​​=0.2cos(9π​n)+2cos(10049π​n)+0.4cos(54π​n)

此时三个频率分量

2 π ω 1 = 2 π × 200 99 π = 400 99 2 π ω 2 = 2 π × 100 49 π = 200 49 2 π ω 1 = 2 π × 5 4 π = 5 2 \begin{array}{l} \frac{{2\pi }}{{{\omega _1}}} = 2\pi \times \frac{{200}}{{99\pi }} = \frac{{400}}{{99}}\\ \frac{{2\pi }}{{{\omega _2}}} = 2\pi \times \frac{{100}}{{49\pi }} = \frac{{200}}{{49}}\\ \frac{{2\pi }}{{{\omega _1}}} = 2\pi \times \frac{5}{{4\pi }} = \frac{5}{2} \end{array} ω1​2π​=2π×99π200​=99400​ω2​2π​=2π×49π100​=49200​ω1​2π​=2π×4π5​=25​​

由此可知最小公倍数为N=400，所以x(n)的周期为N=400。

一个周期的数据长度为 T 0 = N T s = 400 × 0.25 m s = 0.1 s {T_0} = N{T_s} = 400 \times 0.25ms = 0.1s T0​=NTs​=400×0.25ms=0.1s

当抽样点数为周期的整数倍时便不会产生频谱泄露。

由于信号中最小的频率间隔为 ( f 1 − f 2 = 10 H z ) ({f_1} - {f_2} = 10Hz) (f1​−f2​=10Hz)所以 F 0 = 10 H z F_0=10Hz F0​=10Hz，为此所需的数据长度为 T 0 = 1 F 0 = 1 10 H z = 0.1 s {T_0} = \frac{1}{{{F_0}}} = \frac{1}{{10Hz}} = 0.1s T0​=F0​1​=10Hz1​=0.1s

与一个周期的数据长度相同，可取抽样点数 N ≥ f s f 1 − f 2 = 400 N \ge \frac{{{{\rm{f}}_{\rm{s}}}}}{{{f_1} - {f_2}}} = 400 N≥f1​−f2​fs​​=400,因此最少需要的采样点数N取400。

数字角频率 ω = 2 π f f S \omega = \frac{{2\pi f}}{{{f_S}}} ω=fS​2πf​,得数字频率分别为 ω 1 = 0.154 , ω 2 = 0.156 , ω 3 = 0.251 {\omega _1} = 0.154,{\omega _2} = 0.156,{\omega _3} = 0.251 ω1​=0.154,ω2​=0.156,ω3​=0.251。

采用MATLAB方法做频谱分析,结果如图。

clear all;clc;Ts=0.00025;%采样间隔为0.25msfs=1/Ts;N=400; %抽样点N取400，抽样频率fs取4000n=[0:1:N-1];k=[0:1:N-1];xn=0.2*cos(99*(pi/200)*n)+2*cos(49*(pi/100)*n)+0.4*cos(4*(pi/5)*n);%信号离散序列WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.^nk;Xk=xn*WNnk; %根据定义求信号DFTX=abs(Xk);figure(1);stem(k,X,'*');grid on;xlabel('k');title('幅频特性曲线（K=400，fs=4000Hz）');ylabel('|X(k)|');figure(2);stem(2*pi/fs*k,X,'*');grid on;xlabel('数字频率');title('数字频率（K=400，fs=4000Hz）');ylabel('|X(k)|');figure(3);stem(fs/N*k,X,'*');grid on;xlabel('模拟频率');title('模拟频率（K=400，fs=4000Hz）');ylabel('|X(k)|');

频谱泄漏

（3）如果采样点数N=512，样本信号 x ( n T s ) x(nT_s) x(nTs​)的DFT中，是否存在频谱泄漏，如存在，请分析造成频谱泄露的原因、造成的可能影响，以及可能缓轻它的方法？

解：分析可知当采样点数N为512时，此时的频谱 k = 125 时 , f 1 = k N × f s = 125 512 × 4000 = 976.6 H z ; k = 126 时 , f 1 = k N × f s = 126 512 × 4000 = 984.4 H z k = 125时,{f_1} = \frac{k}{N} \times {f_s} = \frac{{125}}{{512}} \times 4000 = 976.6Hz;k = 126时,{f_1} = \frac{k}{N} \times {f_s} = \frac{{126}}{{512}} \times 4000 = 984.4Hz k=125时,f1​=Nk​×fs​=512125​×4000=976.6Hz;k=126时,f1​=Nk​×fs​=512126​×4000=984.4Hz此时980Hz的频率分量在此两频率之间。

k = 126 时， f 1 = k N × f s = 126 512 × 4000 = 984.4 H z ; k = 127 时， f 1 = k N × f s = 127 512 × 4000 = 992.2 H z k = 126时，{f_1} = \frac{k}{N} \times {f_s} = \frac{{126}}{{512}} \times 4000 = 984.4Hz;k = 127时，{f_1} = \frac{k}{N} \times {f_s} = \frac{{127}}{{512}} \times 4000 = 992.2Hz k=126时，f1​=Nk​×fs​=512126​×4000=984.4Hz;k=127时，f1​=Nk​×fs​=512127​×4000=992.2Hz此时990Hz的频率分量在此两频率之间。

k = 204 时， f 1 = k N × f s = 204 512 × 4000 = 1593.8 H z ; k = 205 时， f 1 = k N × f s = 205 512 × 4000 = 1601.6 H z k = 204时，{f_1} = \frac{k}{N} \times {f_s} = \frac{{204}}{{512}} \times 4000 = 1593.8Hz;k = 205时，{f_1} = \frac{k}{N} \times {f_s} = \frac{{205}}{{512}} \times 4000 = 1601.6Hz k=204时，f1​=Nk​×fs​=512204​×4000=1593.8Hz;k=205时，f1​=Nk​×fs​=512205​×4000=1601.6Hz此时1600Hz的频率分量在此两频率之间。

此时的频率分辨率 F 0 = f s N = 7.8 < 10 H z {F_0} = \frac{{{f_s}}}{N} = 7.8 < 10Hz F0​=Nfs​​=7.8<10Hz能够满足分辨率的要求，但是会出现少量的频谱泄漏。

造成频谱泄漏的原因：所取的序列点数不是序列周期的整数倍。频谱泄漏的产生与信号时域截断的长度有关，因为在数字信号处理中，信号长度不可能无限长，信号长度越长，分辨率越高，频谱泄漏越少。因为对时域信号进行截断，本质上就是乘了一个矩形窗，也就是相当于在频率域做了卷积。

可能造成的影响：频谱泄露的直接影响是无法准确获得频率幅值，因为在某频率的信号幅值被周围的频率值平分。截短后信号不能准确代表原始信号，谱分析产生误差，出现了不该有的频率分量，从而导致对于频谱的分析不准确。产生的较多的旁瓣还有可能带来谱间干扰和混叠失真。原有谱线被展宽，增加了不必要的频率分量。降低频率分辨率（矩形窗谱主瓣宽度的一半）。谱间干扰，强信号的旁瓣影响弱信号的主瓣（淹没信号）。

缓轻频谱泄漏方法：

1、对周期信号——采样和分析的点数为周期的整倍数；

2、减小主瓣宽度，主瓣宽度与截取长度N成反比，即增加x(n)长度；

3、减小旁瓣宽度，加合适的窗函数，缓慢截断；

4、整周期截断——截取的正好是序列的一个或整数个周期。可以代表原序列，不发生频谱泄漏。

使用MATLAB作采样点N=512时的频谱分析，如图发生了频谱泄漏。

clear all;clc;Ts=0.00025;%采样间隔为0.25msfs=1/Ts;N=512; %抽样点N取512，抽样频率fs取4000n=[0:1:N-1];k=[0:1:N-1];xn=0.2*cos(99*(pi/200)*n)+2*cos(49*(pi/100)*n)+0.4*cos(4*(pi/5)*n);%信号离散序列WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.^nk;Xk=xn*WNnk; %根据定义求信号DFTX=abs(Xk);subplot(2,1,1);stem(k,X,'*');grid on;xlabel('k');title('幅频特性曲线（K=512，fs=4000Hz）');ylabel('|X(k)|');x=fs/N*k;subplot(2,1,2);stem(x,X,'*');grid on;xlabel('模拟频率');title('模拟频率（K=512，fs=4000Hz）');ylabel('|X(k)|');

FIR数字滤波器设计

（4）如希望在样本 x ( n T s ) x(nT_s) x(nTs​)中，通过FIR数字滤波器有效的抑制 x 3 ( t ) x_3(t) x3​(t)的频谱分量，使其衰减超过50dB。试设计该数字滤波器（方法不限），给出设计指标及设计方案，并对所设计的滤波器进行仿真，分析其性能。

解：由题意可知，信号 x 1 ( t ) x_1(t) x1​(t)、 x 2 ( t ) x_2(t) x2​(t)和 x 3 ( t ) x_3(t) x3​(t)的频率分别为0.99kHz，0.98kHz和1.6kHz，为有效的抑制 x 3 ( t ) x_3(t) x3​(t)的频谱分量，使其衰减超过50dB。结合信号频率情况此处令 f p = 1.1 k H z , f s t = 1.5 k H z , A s = 50 d B {f_p} = 1.1kHz,{f_{st}} = 1.5kHz,{A_s} = 50dB fp​=1.1kHz,fst​=1.5kHz,As​=50dB。

考虑采用窗函数设计法。

窗函数法设计线性相位FIR DF的设计步骤为:

1、确定所需滤波器的性能指标要求。

2、按4种理想线性相位DF(低通、带通、带阻,高通)的表达式,将理想滤波器的截止频率(或上下截止频率)确定为过渡带的算术中心频率。

3、得到理想滤波器的单位冲激响应 h d ( n ) h_d(n) hd​(n)。

4、选择窗函数的类型和长度点数N。窗函数的类型由给定滤波器阻带最小衰减As(dB)确定。窗的长度点数N由滤波器过渡带宽来确定。

5、求加窗后的实际滤波器的单位冲激响应 h ( n ) = h d ( n ) ω ( n ) h(n) = {h_d}(n)\omega (n) h(n)=hd​(n)ω(n)。

6、检验 H ( e j ω ) = D T F T [ h ( n ) ] H({e^{j\omega }}) = DTFT[h(n)] H(ejω)=DTFT[h(n)]是否满足所求滤波器的性能要求。

将模拟频率转换为数字频率： ω p = 2 π f p / f s = 0.55 π ω s t = 2 π f s t / f s = 0.75 π \begin{array}{l} {\omega _p} = 2\pi {f_p}/{f_s} = 0.55\pi \\ {\omega _{st}} = 2\pi {f_{st}}/{f_s} = 0.75\pi \end{array} ωp​=2πfp​/fs​=0.55πωst​=2πfst​/fs​=0.75π​

理想低通滤波波器的频率响应 H d ( e j ω ) {H_d}\left( {{e^{j\omega }}} \right) Hd​(ejω)和单位抽样响应 h d ( n ) {h_d}(n) hd​(n)为：

其中 τ = ( N − 1 ) / 2 \tau = (N - 1)/2 τ=(N−1)/2， ω c = 0.5 ω p + ω s t = 0.65 π {\omega _c} = 0.5{\omega _p} + {\omega _{st}} = 0.65\pi ωc​=0.5ωp​+ωst​=0.65π。

由要求的阻带衰减要超过 A s = 50 d B {A_s} = 50dB As​=50dB，由于海明窗阻带最小衰减为53dB，故选择窗函数的类型为海明窗。

过渡带宽 Δ ω = ω s t − ω p = 0.75 π − 0.55 π = 0.2 π \Delta \omega = {\omega _{st}} - {\omega _p} = 0.75\pi - 0.55\pi = 0.2\pi Δω=ωst​−ωp​=0.75π−0.55π=0.2π。通过查表，海明窗过渡带为 Δ ω = 6.6 π / N \Delta \omega = 6.6\pi /N Δω=6.6π/N

所以可知 N = 6.6 π Δ ω = 6.6 π 0.2 π = 33 N = \frac{{6.6\pi }}{{\Delta \omega }} = \frac{{6.6\pi }}{{0.2\pi }} = 33 N=Δω6.6π​=0.2π6.6π​=33

由此群延迟 τ \tau τ： τ = N − 1 2 = 16 \tau = \frac{{N - 1}}{2} = 16 τ=2N−1​=16

因此，海明窗 ω ( n ) \omega (n) ω(n)为： ω ( n ) = [ 0.54 − 0.46 cos ⁡ ( 2 π n N − 1 ) ] R N ( n ) = [ 0.54 − 0.46 cos ⁡ π n 16 ] R 33 ( n ) \omega (n) = \left[ {0.54 - 0.46\cos (\frac{{2\pi n}}{{N - 1}})} \right]{R_N}(n){\rm{ = }}\left[ {0.54 - 0.46\cos \frac{{\pi n}}{{16}}} \right]{R_{33}}(n) ω(n)=[0.54−0.46cos(N−12πn​)]RN​(n)=[0.54−0.46cos16πn​]R33​(n)

则线性相位FIR低通滤波器的h(n)为

接下来验证 H ( e j ω ) = Σ n = 0 32 h ( n ) e − j ω n H\left( {{e^{j\omega }}} \right) = \mathop \Sigma \limits_{n = 0}^{32} h(n){e^{ - j\omega n}} H(ejω)=n=0Σ32​h(n)e−jωn。此处采用MATLAB工具画出滤波器幅度响应如图。

由此可以看出阻带衰减还没有达到50dB的要求。所以考虑改选N=35的海明窗进行验证。

ω ( n ) = [ 0.54 − 0.46 cos ⁡ ( 2 π n N − 1 ) ] R N ( n ) = [ 0.54 − 0.46 cos ⁡ π n 17 ] R 35 ( n ) \omega (n) = \left[ {0.54 - 0.46\cos (\frac{{2\pi n}}{{N - 1}})} \right]{R_N}(n){\rm{ = }}\left[ {0.54 - 0.46\cos \frac{{\pi n}}{{17}}} \right]{R_{35}}(n) ω(n)=[0.54−0.46cos(N−12πn​)]RN​(n)=[0.54−0.46cos17πn​]R35​(n)

绘制h(n)和滤波器幅度响应、相位响应如图。

根据程序中对阻带最小衰减的判断确定，As最小为53.22dB,满足了衰减超过50dB的要求。能够有效的抑制 x 3 ( t ) {x_3}(t) x3​(t)的频谱分量。

clc;clear all;fp=1100;fst=1500;%设置通带和阻带频率Ts=0.00025;%采样间隔为0.25msFs=1/Ts; %采样频率wp=2*pi*fp/Fs;ws=2*pi*fst/Fs;%数字角频率deltaw=ws-wp; %求过度带宽N0=ceil(6.6*pi/deltaw);%计算海明窗长N=N0+mod(N0+1,2);%取N为奇数n=[0:N-1];wd=(hamming(N))';%求海明函数wc=(ws+wp)/2;%理想低通滤波器截至频率hd=ideallp(wc,N);%理想低通的单位冲击响应，调用ideallp将其放于相同路径下h=hd.*wd;[db,mag,pha,grd,w] = freqz_m( h,[1]);%检查设计结果的各项指标（“dB”是负值），调用freqz_m将其放于相同路径下dw=2*pi/1000;Rp=-min(db(1:wp/dw +1));%检查通带最大衰减是否满足要求As=-max(db(ws/dw+1:501));%检查阻带最小衰减是否满足要求figure(1); %绘制海明窗plot(n,wd,'linewidth',2);grid;title('海明窗'); xlabel( 'n'); ylabel( 'w(n) ');axis([0,N,0,1]);figure(2);%绘制幅度响应曲线plot(w/pi,db,'linewidth',2)title('幅度响应(db) ' );xlabel('\omega/\pi');ylabel('20log|H(e^j^\omega)|(dB)');axis([0,1,-80,5]);grid;set(gca,'xtickmode','manual','xtick',[ 0,0.55,0.75,1.0]);%横坐标设置，观察感兴趣的位置，此处选取通带和阻带位置figure(3);%绘制相位响应曲线plot(w/pi,pha,'linewidth',2);axis([0,1,-4,4]);gridtitle('相位响应');xlabel ( '\omega/\pi'); ylabel ( 'arg20log[H(e^j^\omega)]');

调用的ideallp和freqz_m函数如下(将他们放于同一工作目录下即可调用)：

function hd = ideallp(wc,N);% 理想低通滤波器的脉冲响应子程序% hd = ideallp(wc,N)% hd = 点0 到 N-1之间的理想脉冲响应% wc = 截止频率(弧度)% N = 理想滤波器的长度%tao = (N-1)/2; % 理想脉冲响应的对称中心位置n = [0:(N-1)]; % 设定脉冲响应长度m = n-tao+eps;% 加一个小数以避免零作除数hd =sin(wc*m) ./ (pi*m); % 理想脉冲响应

function [db,mag,pha,grd,w]=freqz_m(b,a)%求取系统的绝对幅度响应、相对的db值幅度响应、相位响应和群延时响应的函数%db为相对振幅(dB)%mag为绝对振幅%pha为相位响应%grd为群延时%w为频率样本点向量%b为Ha(z)分子多项式系数(对FIR而言，b=h)%a为Hz(z)分母多项式系数(对FIR而言，a=1)% [H,w]=freqz(b,a,1000,'whole');%freqz显示数字滤波器频域中的图形%[H,W] = FREQZ(B,A,N,'whole') uses N points around the whole unit circle.[H,w] = freqz(b,a,1000,'whole');H = (H(1:1:501))'; w = (w(1:1:501))';mag = abs(H);db = 20*log10((mag+eps)/max(mag));pha = angle(H);% pha = unwrap(angle(H));grd = grpdelay(b,a,w);

参考文章

数字信号处理教程第四版_程佩青P210(例3.20) P545 P560(例8.1 例8.7)

链接: 数字信号处理及MATLAB实现

链接: Matlab快速入门（五）傅里叶变换