数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

2024-04-10 1260阅读

温馨提示：这篇文章已超过367天没有更新，请注意相关的内容是否还可用！

数字信号处理翻转课堂笔记18

The Flipped Classroom18 of DSP

对应教材：《数字信号处理（第五版）》西安电子科技大学出版社，丁玉美、高西全著

一、要点

（1）频率采样法设计FIR线性相位滤波器的原理；

（2）线性相位条件对频率响应的约束；

（3）频率采样法设计FIR线性相位滤波器的步骤（重点）；

（4）逼近误差产生的原因及其改进措施（难点，重点）；

（5）基于MATLAB和频率采样法设计FIR线性相位滤波器。

二、问题与解答

1、简述频率采样法设计线性相位FIR滤波器的基本原理。与窗函数法相比，频率采样法具有哪些优势？

2、为什么线性相位条件会制约频率采样法设计FIR滤波器时的频率响应特性？这种制约是否会影响频率采样法所能设计的滤波器类型？

3、自定滤波器长度N、滤波器类型及其理想频率响应特性，按照线性相位条件确定滤波器在[0,2π]区间的频率特性，并在该区间进行N点频域采样，获得频域采样序列H(k)，需注意H(k)应该包含幅度和相位两个部分。然后求H(k)的IDFT，得到对应的滤波器单位响应h(n)。用MATLAB画出h(n)的频率响应（含幅频和相频）波形（连续曲线），验证其是否具有线性相位特性，并与自定的滤波器理想频率特性进行比较，对比其差异。根据频域采样的相关原理，分析差异（即逼近误差）产生的原因。

4、对上一题的频域采样点，在适当的位置加入过渡点，重新绘制幅频特性波形，与上一题的幅频特性进行比较，对比其截止频率、阻带衰减等方面的差异。根据本题和上一题的结果，讨论为什么用频率采样法设计FIR滤波器，通常需要在设计完之后验证其是否满足预期的技术指标要求。

5、基于MATLAB，用频率采样法分别设计一个低通和一个高通滤波器，技术指标自定。在过渡带加入1个或2个过渡点，调整过渡点的幅值大小，画出并比较其幅频特性波形、单位响应波形。通过反复尝试，使得滤波器阻带衰减达到最低，此时过渡点的幅度值是多少？总结频率采样法设计FIR滤波器的步骤，讨论过渡点数量和幅值对滤波器频率特性（包括过渡带宽度、阻带衰减、通带纹波、截止频率等）的影响，对比低通、高通滤波器单位响应的显著差异。（滤波器阶次N怎么确定）

6、利用FIR2函数设计上一题的FIR滤波器（可以只选低通或者高通其中一种），过渡点数量和幅值自定。分别采用3种以上不同的窗函数，画出不同窗函数所设计滤波器的幅频特性，比较其通带纹波、阻带衰减、过渡带宽度、截止频率的差异。讨论采用通过加窗改善幅频特性和通过过渡点改善幅频特性这两种方法各自的特点。

7、x=[-4, -2, 0, -4, -6, -4, -2, -4, -6, -6, -4, -4, -6, -6, -2, 6, 12, 8, 0, -16, -38, -60, -84, -90, -66, -32, -4, -2, -4, 8, 12, 12, 10, 6, 6, 6, 4, 0, 0, 0, 0, 0, -2, -4, 0, 0, 0, -2, -2, 0, 0, -2, -2, -2, -2, 0 ]，这是一个心电图信号，混进了高频噪声，需要设计一个数字低通滤波器，滤除高频噪声，但应保持心电图的波形不失真，请设计滤波器完成对该信号的滤波，比较滤波前后的波形（可参考学堂在线视频6-12，但应修改FIR滤波器设计方法，不要采用等波纹最佳逼近法）。提示：先要对信号进行频谱分析，分析其高频噪声所处的频带，以确定滤波器的截止频率。设计好滤波器之后，对信号的滤波，采用DFT来实现（用循环卷积计算线性卷积），或者直接用conv函数求卷积。如果采用filter函数，应先对x补零一定的长度再进行滤波，否则由于输出序列长度的限制（输出序列长度与x长度一致），可能无法看到滤波后输出的完整波形。

1、频率采样法基本原理

简述频率采样法设计线性相位FIR滤波器的基本原理。与窗函数法相比，频率采样法具有哪些优势？

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

滤波器技术指标通常在频域给出，而窗函数法从时域出发进行滤波器设计，需把频率特性转换为单位脉冲响应。当频率特性表达式较为复杂，或无闭合形式时，单位脉冲响应求取较为困难。

频率采样法直接从频域进行设计，物理概念清楚，直观方便；尤其适合于窄带滤波器，这时频率响应只有少数几个非零值，用频率采样网络架构来实现非常方便。

2、线性相位条件

为什么线性相位条件会制约频率采样法设计FIR滤波器时的频率响应特性？这种制约是否会影响频率采样法所能设计的滤波器类型？

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

会影响，对于高通和带阻滤波器，n只能取奇数。

3、频率采样与线性相位

自定滤波器长度N、滤波器类型及其理想频率响应特性，按照线性相位条件确定滤波器在[0,2π]区间的频率特性，并在该区间进行N点频域采样，获得频域采样序列H(k)，需注意H(k)应该包含幅度和相位两个部分。然后求H(k)的IDFT，得到对应的滤波器单位响应h(n)。用MATLAB画出h(n)的频率响应（含幅频和相频）波形（连续曲线），验证其是否具有线性相位特性，并与自定的滤波器理想频率特性进行比较，对比其差异。根据频域采样的相关原理，分析差异（即逼近误差）产生的原因。

figure(1)
Bt=pi/5;wp=pi/2;
N=27;                   %使N为奇数
Np=fix(wp/(2*pi/N));    %Np+1为通带[0,wp]上采样点数
Ns=N-2*Np-1;            %Ns为阻带[wp,2*pi一wp]上采样点数
Hk=[ones(1,Np+1),zeros(1,Ns),ones(1,Np)];   %N为奇数，设定理想幅度采样向量偶对称
thetak=-pi*(N-1)*(0:N-1)/N;                 %相位采样向量0(k)=N一1)πk/,0≤k≤N一1
Hdk=Hk.*exp(j*thetak);                      %构造频域采样向量Hd(k)
hn=real(ifft(Hdk));                         %h (n)=IDFT [H(k)]
Hw=fft (hn,1024)                            %计算频率响应函数：DFT[h(n)]
wk=2*pi*(0:1023)/1024;
Hgw=Hw.*exp(j*wk*(N-1)/2);                  %计算幅度响应函数计算通带最大纹波RD和阻带最小衰减Rs
Rp=max(20*log10(abs(Hgw)));
hgmin=min(real(Hgw));
Rs=20*log10(abs (hgmin));
stem(hn)
figure(2)
subplot(2,1,1)
plot (wk/pi,abs (Hgw))
hold on
plot(0:1/800:0.25+1/800,[ones(1,201),0],2-0.25-1/800:1/800:2,[0,ones(1,201)]);
hold off
subplot (2,1,2)
plot (wk/pi,angle(Hw))

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

4、加入过渡点的频域采样

对上一题的频域采样点，在适当的位置加入过渡点，重新绘制幅频特性波形，与上一题的幅频特性进行比较，对比其截止频率、阻带衰减等方面的差异。根据本题和上一题的结果，讨论为什么用频率采样法设计FIR滤波器，通常需要在设计完之后验证其是否满足预期的技术指标要求。

T=0.5;                  %输入过渡带采样值T
Bt=pi/5;wp=pi/2;
N=27;
Np=fix(wp/(2*pi/N));    %Np+1为通带[0,wp]上采样点数
Ns=N-2*Np-1;            %Ns为阻带[wp,2*pi一wp]上采样点数
Hk=[ones(1,Np+1),zeros(1,Ns),ones(1,Np)];   %N为奇数，设定理想幅度采样向量偶对称
Hk(Np+2)=T;Hk(N-Np)=T;  %加一个过渡采样
thetak=-pi*(N-1)*(0:N-1)/N;                 %相位采样向量θ(k)=(N-1)πk/N,0≤k≤N-1
Hdk=Hk.*exp(j*thetak);                      %构造频域采样向量Hd(k)
hn=real(ifft(Hdk));                         %h(n)=IDFT[H(k)]
Hw=fft(hn,1024);                            %计算频率响应函数：DFT[h(n)]
wk=2*pi*(0:1023)/1024;
Hgw=Hw.*exp(j*wk*(N-1)/2);                 %计算幅度响应函数计算通带最大纹波R即和阻带最小衰减Rs
Rp=max (20*log10(abs (Hgw)));
hgmin=min(real (Hgw));
Rs=20*log10(abs(hgmin));
figure(3);
stem(hn);
figure(4);
subplot (2,1,1);
plot (wk/pi,abs (Hgw));
hold on;
plot(0:1/800:0.25+1/800,[ones(1,201),0],2-0.25-1/800:1/800:2,[0,ones(1,201)]);
hold off
subplot (2,1,2)
plot (wk/pi,angle(Hw))

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

讨论：截止频率变大，阻带衰减变大,直接对理想滤波器的频率响应采样的“基本频率采样法”不能满足一般工程对阻带衰减的要求。

5、加入过渡点对滤波器频率特性的影响

基于MATLAB，用频率采样法分别设计一个低通和一个高通滤波器，技术指标自定。在过渡带加入1个或2个过渡点，调整过渡点的幅值大小，画出并比较其幅频特性波形、单位响应波形。通过反复尝试，使得滤波器阻带衰减达到最低，此时过渡点的幅度值是多少？总结频率采样法设计FIR滤波器的步骤，讨论过渡点数量和幅值对滤波器频率特性（包括过渡带宽度、阻带衰减、通带纹波、截止频率等）的影响，对比低通、高通滤波器单位响应的显著差异。（滤波器阶次N怎么确定）

T=0.2;                     
%不断更改T值
Bt=pi/18;wp=pi/3;
m=1;
N=ceil((m+1)*2*pi/Bt);      %估算采样点数N
N=N+mod(N+1,2);             %使N为奇数
Np=fix(wp/(2*pi/N));        %Np+1为通带[0,wp]上的采样点数
Ns=N-2*Np-1;                %Ns为阻带[wp,2*pi-wp]上的采样点数
Hk=[ones(1,Np+1),zeros(1,Ns),ones(1,Np)];   %N为奇数，幅度采样采样向量偶对称A(k)-A(N-k)
Hk(Np+2)=T; Hk(N-Np)=T;                     %加一个过渡采样
thetak=-pi*(N-1)*(0:N-1)/N;     %相位采样向量
Hdk=Hk.*exp(j*thetak);          %构造频域采样向量Hd(k)
hn=real(ifft(Hdk));
N=length(hn);
Hw=fft(hn,1024);                %计算频率响应函数
wk=2*pi*[0:1023]/1024;
Hgw=Hw.*exp(j*wk*(N-1)/2);      %计算幅度响应函数
Rp=max(20*log10(abs(Hgw)));     %计算通带最大衰减RP
hgmin=min(real(Hgw));
Rs=20*log10(abs(hgmin));        %计算阻带最小衰减
subplot(2,1,1);
plot(wk/pi,20*log10(abs(Hgw)));
xlabel('w/pi');
ylabel('20*log10(abs(Hgw)))');
xlabel('1个采样点，T=0.2');
subplot(2,1,2);
n=0:N-1;
stem(n,hn);
xlabel('n');
ylabel('hn');

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

T=0.38

Rp = 0.4797 Rs = -43.3937

T=0.5

Rp = 0.2802 Rs = -29.7060

T=0.2

Rp = 0.7706 Rs = -30.2164

高通：

T=0.2;
%不断更改T
Bt=pi/18;wp=pi*0.7;
m=1;
N=ceil((m+1)*2*pi/Bt);%估算采样点数N
N=N+mod(N+1,2);  %使N为奇数
Np=fix(wp/(2*pi/N)); %Np+1为通带[0,wp]上的采样点数
Ns=N-2*Np-1; %Ns为阻带[wp,2*pi-wp]上的采样点数
Hk=[zeros(1,Np+1),ones(1,Ns),zeros(1,Np)]; %N为奇数，幅度采样采样向量偶对称A(k)=A(N-k)
Hk(Np+2)=T; Hk(N-Np)=T;                     %加一个过渡采样
thetak=-pi*(N-1)*(0:N-1)/N;     %相位采样向量
Hdk=Hk.*exp(j*thetak);          %构造频域采样向量Hd(k)
hn=real(ifft(Hdk));
N=length(hn);
Hw=fft(hn,1024);                %计算频率响应函数
wk=2*pi*[0:1023]/1024;
Hgw=Hw.*exp(j*wk*(N-1)/2);      %计算幅度响应函数
Rp=max(20*log10(abs(Hgw)));     %计算通带最大衰减RP
hgmin=min(real(Hgw));
Rs=20*log10(abs(hgmin));        %计算阻带最小衰减
subplot(2,1,1);
plot(wk/pi,20*log10(abs(Hgw)));
xlabel('w/pi');
ylabel('20*log10(abs(Hgw)))');
xlabel('1个采样点，T=0.2');
subplot(2,1,2);
n=0:N-1;
stem(n,hn);
xlabel('n');
ylabel('hn');

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

T=0.35

Rp = 0.5687 Rs = -46.5289

T=0.5

Rp = 0.2823 Rs = -29.6693

T=0.2

Rp = 0.8500 Rs = -33.3512

多个采样点（仅供参考）

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

T1=0.5 T2=0.06

Rp = 0.3350 Rs = -53.0071

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

T1=0.06 T2=0.5

Rp = 0.3537 Rs = -63.3920

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

当过渡点的幅值小于最佳幅值时，截止频率减小，通带最大衰减较小，阻带最小衰减降低，其阻带的纹波下降比较快，当过渡点的幅值大于最佳幅值时，通带最大衰减较大，阻带最小衰减较小，截止频率也小于最佳幅值情况下的截止频率，过渡带宽度较小。

高通滤波器单位响应与低通滤波器的单位响应之间的差异：高通滤波器通带频率较高，反映在波形上为每个点位正负相间，但由于有频率较为低的部分，所以有的相邻点之间均为正或均为负，对于低通滤波器，由于通带频率较低，符号变化没有高通的快。

过渡点数量增多时，滤波器阻带最小衰减增大，但过渡带带宽也增大。（未验证）

6、不同窗函数的滤波特性

利用FIR2函数设计上一题的FIR滤波器（可以只选低通或者高通其中一种），过渡点数量和幅值自定。分别采用3种以上不同的窗函数，画出不同窗函数所设计滤波器的幅频特性，比较其通带纹波、阻带衰减、过渡带宽度、截止频率的差异。讨论采用通过加窗改善幅频特性和通过过渡点改善幅频特性这两种方法各自的特点。

T1=0.38;  %输入过渡带采样值T
Bt=pi/18;wp=pi/3;%过渡带宽度pi/13，通带截止频率pi/4
m=1; %过渡点个数
N=ceil((m+1)*2*pi/Bt)+1;  %估算采样点数N
F=[0,wp/pi,wp/pi+2/N,wp/pi+4/N,1];A=[1,1,T1,0,0];
hn=fir2(N-1,F,A,boxcar(N));
N=length(hn);
n=0:N-1;
Hw=fft(hn,1024); %计算频率响应函数wk=2*pi*[0:1023]/1024;
wk=2*pi*[0:1023]/1024;
Hgw=Hw.*exp(j*wk*(N-1)/2); %计算幅度响应函数
Rp=max(20*log10(abs(Hgw)))%计算通带最大衰减RP
hgmin=min(real(Hgw));
Rs=20*log10(abs(hgmin)) %计算阻带最小衰减
subplot(3,2,1)
stem(n,hn)
xlabel('n')
ylabel('hn')
subplot(3,2,2)
plot(wk/pi,20*log10(abs(Hgw)));
xlabel('w/\pi')
ylabel('20*log10(abs(Hgw)) dB')
title('矩形窗')
grid on
 
hn=fir2(N-1,F,A,hanning(N));
N=length(hn);
n=0:N-1;
Hw=fft(hn,1024); %计算频率响应函数wk=2*pi*[0:1023]/1024;
wk=2*pi*[0:1023]/1024;
Hgw=Hw.*exp(j*wk*(N-1)/2); %计算幅度响应函数
Rp=max(20*log10(abs(Hgw)))%计算通带最大衰减RP
hgmin=min(real(Hgw));
Rs=20*log10(abs(hgmin)) %计算阻带最小衰减
subplot(3,2,3)
stem(n,hn)
xlabel('n')
ylabel('hn')
subplot(3,2,4)
plot(wk/pi,20*log10(abs(Hgw)));
xlabel('w/\pi')
ylabel('20*log10(abs(Hgw)) dB')
title('汉宁窗')
grid on
 
hn=fir2(N-1,F,A,blackman(N));
N=length(hn);
n=0:N-1;
Hw=fft(hn,1024); %计算频率响应函数wk=2*pi*[0:1023]/1024;
wk=2*pi*[0:1023]/1024;
Hgw=Hw.*exp(j*wk*(N-1)/2); %计算幅度响应函数
Rp=max(20*log10(abs(Hgw)))%计算通带最大衰减RP
hgmin=min(real(Hgw));
Rs=20*log10(abs(hgmin)) %计算阻带最小衰减
subplot(3,2,5)
stem(n,hn)
xlabel('n')
ylabel('hn')
subplot(3,2,6)
plot(wk/pi,20*log10(abs(Hgw)));
xlabel('w/\pi')
ylabel('20*log10(abs(Hgw)) dB')
title('布莱克曼窗')
grid on

数字信号处理翻转课堂笔记18——频率采样法设计FIR滤波器及matlab实现

通带衰减：矩形窗

VPS购买请点击我

免责声明：我们致力于保护作者版权，注重分享，被刊用文章因无法核实真实出处，未能及时与作者取得联系，或有版权异议的，请联系管理员，我们会立即处理! 部分文章是来自自研大数据AI进行生成,内容摘自(百度百科,百度知道,头条百科,中国民法典,刑法,牛津词典,新华词典,汉语词典,国家院校,科普平台)等数据,内容仅供学习参考,不准确地方联系删除处理! 图片声明：本站部分配图来自人工智能系统AI生成,觅知网授权图片,PxHere摄影无版权图库和百度，360，搜狗等多加搜索引擎自动关键词搜索配图，如有侵权的图片，请第一时间联系我们，邮箱：ciyunidc@ciyunshuju.com。本站只作为美观性配图使用,无任何非法侵犯第三方意图,一切解释权归图片著作权方,本站不承担任何责任。如有恶意碰瓷者,必当奉陪到底严惩不贷!

相关阅读

目录[+]