FIR数字滤波器设计与软件实现

1、设计内容
（1）调用信号产生函数xtg产生具有加性噪声的信号xt，并自动显示xt及其频谱；
（2）设计低通滤波器，从高频噪声中提取xt中的单频抑制载波调幅信号，要求信号幅频失真小于0.1dB，将噪声频谱衰减60dB。观察xt的频谱，确定滤波器指标参数。
（3）根据滤波器指标选择合适的窗函数，计算窗函数的长度N，调用MATLAB函数firl设计一个FIR低通滤波器，并编写程序，调用MATLAB快速卷积函数fftfilt实现对xt的滤波。绘图显示滤波器的频响特性曲线、滤波器输出信号的幅频特性图和时域波形图。
（4）重复（3），滤波器指标不变，但改用等波纹最佳逼近设计FIR滤波器，调用MATLAB函数remezord和remez设计FIR数字滤波器。比较两种设计方法的滤波器的阶数。

function xt=xtg(N)

%信号产生函数，并显示信号的时域波形和幅频特性曲线
%xt=xtg产生一个长度为N，有加性高频噪声的单频调幅信号xt，N=1000；
%采样频率Fs=1000Hz
%载波频率fc=Fs/10=100Hz，调制正弦波频率f0=fc/10=10Hz.
N=1000;Fs=1000;T=1/Fs;Tp=N*T;
t=0:T:(N-1)*T;
fc=Fs/10;f0=fc/10;  %载波频率fc=Fs/10，单频调制信号频率为f0=Fc/10;
mt=cos(2*pi*f0*t);  %产生单频正弦波调制信号mt，频率为f0
ct=cos(2*pi*fc*t);  %产生载波正弦波信号ct，频率为fc
xt=mt.*ct;          %相乘产生单频调制信号xt
nt=0.2*rand(1,N)-1;   %产生随机噪声nt

%======= 设计高通滤波器hn,用于滤除噪声nt中的低频成分,生成高通噪声 =======
fp=150; fs=200;Rp=0.1;As=70;  %滤波器指标
fb=[fp,fs];m=[0,1]; 		  %计算remezord函数所需参数f,m,dev
dev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)];
[n,fo,mo,W]=remezord(fb,m,dev,Fs);	%确定remez函数所需参数
hn=remez(n,fo,mo,W);	    %调用remez函数进行设计,用于滤除噪声nt中的低频成分
yt=filter(hn,1,10*nt);      %滤除随机噪声中低频成分，生成高通噪声yt

%=========以下为绘图部分=========
xt=xt+yt;           
fst=fft(xt,N);k=0:N-1;f=k/Tp;
subplot(2,1,1);
plot(t,xt);
grid on;
xlabel('t/s');ylabel('x(t)');
axis([0,Tp/5,min(xt),max(xt)]);
title('(a) 信号加噪声波形')
subplot(2,1,2);
plot(f,abs(fst)/max(abs(fst)));
grid on;
xlabel('f/Hz');ylabel('幅度')
axis([0,Fs/2,0,1.2]);
title('(b) 信号加噪声频谱')

clear;
%调用xtg产生具有加性噪声的信号x(t)
%x(t)长度N=1000,并显示x(t)信号波形及其频谱
N=1000;xt=xtg;Fs=1000;T=1/Fs;Tp=N*T;
k=0:N-1;f1=k/Tp; 
%初始化指标
%通带截止频率fp 阻带截止频率fs 通带最大衰减 阻带最小衰减 采样频率Fs=1000Hz
fp=120;fs=150;Rp=0.2;As=60;Fs=1000;  

%用窗函数法设计滤波器
wc=(fp+fs)/Fs;      %理想低通滤波器截止频率(关于pi归一化）
B=2*pi*(fs-fp)/Fs;  %过渡带宽度指标
Nw=ceil(11*pi/B);   %blackman窗的长度N ceil取最大整数
hn=fir1(Nw-1,wc,blackman(Nw)); %阻带最小衰减为60dB 可以用blackman窗
disp(Nw-1);             %窗函数法设计滤波器 滤波器长度-1=阶数
Hw=abs(fft(hn,1024));   %求设计的滤波器频率特性
ywt=fftfilt(hn,xt,N);   %调用函数fftfilt对xt滤波

%窗函数法设计法的绘图部分（滤波器损耗函数, 滤波器输出信号波形）
f=(0:1023)*Fs/1024;
figure(2);subplot(3,1,1);
plot(f,20*log10(Hw/max(Hw)),'g');
grid on;
title('窗函数法低通滤波器幅频特性')
axis([0,Fs/2,-120,20]);
xlabel('f/Hz');ylabel('幅度')
t=(0:N-1)/Fs;Tp=N/Fs; 
subplot(3,1,2);
plot(t,ywt,'r');grid;
axis([0,Tp/2,-1,1]);xlabel('t/s');ylabel('y_w(t)');
title('滤除噪声后的信号波形');

subplot(3,1,3);
fxt=fftshift(fft(ywt,1000));
stem(f1-500,abs(fxt)/max(abs(fxt)),'k.');
grid on;
title('滤除噪声后的信号频谱');
axis([0,500,0,1]);      %设置刻度  右边

%等波纹最佳逼近法设计滤波器
fb=[fp,fs]; m=[1,0];        
%确定remezord函数所需参数f,m,dev
dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)]; %十进制表示的
[Ne,fo,mo,W]=remezord(fb,m,dev,Fs);
disp(Ne-1);   %波纹最佳逼近法设计滤波器 滤波器长度
%确定remez函数所需参数
hn=remez(Ne-1,fo,mo,W);  %调用remez函数进行设计
Hw=abs(fft(hn,1024));  %求设计的滤波器频率特性
yet=fftfilt(hn,xt,N);  %调用函数fftfilt对xt滤波

%等波纹设计法的绘图部分（滤波器损耗函数, 滤波器输出信号yw(t)波形）
f=(0:1023)*Fs/1024;
figure(3);subplot(3,1,1);
plot(f,20*log10(Hw/max(Hw)),'g');
grid on;
title('等波纹法低通滤波器幅频特性');
axis([0,Fs/2,-80,10]);
xlabel('f/Hz');ylabel('幅度')
t=(0:N-1)/Fs;Tp=N/Fs;
subplot(3,1,2)
plot(t,yet,'r');grid;
axis([0,Tp/2,-1,1]);xlabel('t/s');ylabel('y_e(t)');
title('滤除噪声后的信号波形');

subplot(3,1,3);
fet=fftshift(fft(yet,1000));
stem(f1-500,abs(fet)/max(abs(fet)),'k.');title('滤除噪声后的信号频谱');
axis([0,500,0,1]);      %设置刻度

2、思考题
（1）如果要求用窗函数法设计带通滤波器，给定通带上、下截止频率和阻带上、下截止频率，如何确定理想带通滤波器的截止频率。
wc=(fp+fs)/Fs;
（2）解释为什么对同样的技术指标，用等波纹最佳逼近法设计的滤波器阶数低。
因为窗函数法通带最大衰减阻带最大衰减不能控制，每一个窗都有固定的最大衰减，存在较大富裕，等波纹法，通带和阻带均为等波纹特性，最大衰减可控，没有资源浪费，所以在相同指标下阶数较低。