如果你只是要理论仿真的话,切比雪夫逼近法最好,FDATOOL就行了,如果只想满足自己想法用频率抽样法在Matlab里溜一遍,我做过用频率采样发设计的FIR低通滤波器,高通自己改下参数和符号数,程序说明很清楚了,参考胡广书的《数字信号处理-理论、算法与实现》改程序里的符号变量和数据。
代码:
%-----------------------------------------------------------------
%作者:樊彦(10信处)
%名称:数字信号处理---胡广书
%功能:习题7.3——FIR设计数字低通滤波器(DLP)→EvenOdd
%理想低通频率响应:Hd=1:|w|<=pi/4;Hd=0:pi/4<|w|<pi;
%版本:3.0
%-----------------------------------------------------------------
clear all;
clc;
N=51; %滤波器系数的长度(时域点数)
M=N-1; %辅助参数
Wc=pi/4; %理想截止频率
%-----------------------------------------------------------------
%符号运算,求hn→FIR-DLP的时域序列
%--------------------------------
syms k n real %指定符号变量→w:频率,n:离散时间
%N为偶数时
Hd_Even_Front=exp(-j*k*pi*(N-1)/N); %0,1,...,N/2-1
Hd_Even_Half=0; %N/2
Hd_Even_Back=-exp(-j*k*pi*(N-1)/N); %N/2+1,N/2+2,...,N-1
%N为奇数时
Hd_Odd=exp(-j*k*pi*(N-1)/N); %0,1,...,N-1
%-----------------------------------------------------------------
%先计算Hd,再计算hn→IDFT(Hd)
nc=floor(N*Wc/(2*pi)); %计算通带Hd截止(Matlab)坐标
Hd_value=0.3830; %改进滤波器性能参数(0.3~0.5)
if 0~=rem(N,2) %N为奇数
for i=0:N-1
if i<=nc
Hd(i+1)=double(limit(Hd_Odd,k,i)); %计算Hd_奇数
elseif i>=N-nc
Hd(i+1)=double(limit(Hd_Odd,k,i)); %计算Hd_奇数
else
Hd(i+1)=0;
end
end
Hd(1,nc+2)=Hd_value*double(limit(Hd_Odd,k,nc+1)); %改善Hd—1
Hd(1,N-nc)=Hd_value*double(limit(Hd_Odd,k,N-nc-1));
% Hd(1,nc+3)=0.01*double(limit(Hd_Odd,k,nc+2)); %改善Hd—2
% Hd(1,N-nc-1)=0.01*double(limit(Hd_Odd,k,N-nc-2));
else %N为偶数
for i=0:N-1
if i<=(N/2-1)
if i<=nc
Hd(i+1)=double(limit(Hd_Even_Front,k,i)); %计算Hd_偶数
else
Hd(i+1)=0;
end
elseif i==N/2
Hd(i+1)=0; %计算Hd_偶数
else
if i>=N-nc
Hd(i+1)=double(limit(Hd_Even_Back,k,i)); %计算Hd_偶数
else
Hd(i+1)=0;
end
end
end
Hd(1,nc+2)=Hd_value*double(limit(Hd_Even_Front,k,nc+1)); %改善Hd—1
Hd(1,N-nc)=Hd_value*double(limit(Hd_Even_Back,k,N-nc-1));
% Hd(1,nc+3)=0.01*double(limit(Hd_Even_Front,k,nc+2)); %改善Hd—2
% Hd(1,N-nc-1)=0.01*double(limit(Hd_Even_Back,k,N-nc-2));
end
hn=ifft(Hd,N); %计算hn
%-----------------------------------------------------------------
%求频率响应
%-------------------
%b=hn; %FIR滤波器的系数
Hw_0=sum(hn); %w=0处的IDTFT(hn)
b=hn/Hw_0; %FIR滤波器的归一化系数
a=1;
[H,w]=freqz(b,a); %返回w为数字频率值
%-----------------------------------------------------------------
%幅频响应&&相频响应
Hr=abs(H); %幅频响应
dBHr=20*log10(Hr); %幅频响应(dB)
%dBHr=Hr;
Hphase=angle(H); %相频响应
Hphase=unwrap(Hphase); % 解卷绕
%-----------------------------------------------------------------
%图形-数据可视化
set(figure(1),'Name','时域实数序列图')
stem(0:M,hn,'.');grid on;
axis([-0.1 M+0.1 -1.1*max(abs(hn)) 1.1*max(abs(hn))]); %限定横坐标和纵坐标的显示范围
title('FIR-频率抽样法-DLPF-时域实数序列');
%semilogy(w,Hr);grid on; %x轴线性,y轴对数
ylabel(' Amplitude value ');
xlabel(' Discrete Time (s)');
%-----------------------------------------------------------------
set(figure(2),'Name','频域序列序列图')
subplot(311)
stem(0:M,real(Hd),'.');grid on; %画频率响应的实部序列--关于中点偶对称
axis([-0.1 M+0.1 -1.1*max(abs(real(Hd))) 1.1*max(abs(real(Hd)))]); %限定横坐标和纵坐标的显示范围
title('FIR-频率抽样法-DLPF-频域值的实部序列(关于中点偶对称)');
%semilogy(w,Hr);grid on; %x轴线性,y轴对数
ylabel(' Amplitude value ');
xlabel(' Discrete Time (s)');
subplot(312)
stem(0:M,imag(Hd),'.');grid on; %画频率响应的虚部序列--关于中点奇对称
axis([-0.1 M+0.1 -1.1*max(abs(imag(Hd))) 1.1*max(abs(imag(Hd)))]); %限定横坐标和纵坐标的显示范围
title('FIR-频率抽样法-DLPF-频域值的虚部序列(关于中点奇对称)');
%semilogy(w,Hr);grid on; %x轴线性,y轴对数
ylabel(' Amplitude value ');
xlabel(' Discrete Time (s)');
subplot(313)
stem(0:M,abs(Hd),'.');grid on; %画频率响应的幅值序列--关于中点偶对称
axis([-0.1 M+0.1 -1.1*max(abs(Hd)) 1.1*max(abs(Hd))]); %限定横坐标和纵坐标的显示范围
title('FIR-频率抽样法-DLPF-频域摸值序列');
%semilogy(w,Hr);grid on; %x轴线性,y轴对数
ylabel(' Amplitude value ');
xlabel(' Discrete Time (s)');
%-----------------------------------------------------------------
set(figure(3),'Name','相频响应图')
plot(w,Hphase,'');grid on; %画相频响应虚线
title('FIR-频率抽样法-DLPF-相频响应');
ylabel(' Phase Freq. Res.');
xlabel(' Digital Freq (rad)');
%-----------------------------------------------------------------
set(figure(4),'Name','幅频响应(dB)图')
%subplot(211)
plot(w,dBHr,'');grid on; %画归一化的幅频率响应曲线
axis([0 pi 1.1*min(dBHr) 1.01*max(dBHr)]); %限定横坐标和纵坐标的显示范围
%axis([0 pi/4 -5 1.01*max(dBHr)]);
title('FIR-频率抽样法-DLPF-幅频响应(dB)');
%semilogy(w,Hr);grid on; %x轴线性,y轴对数
ylabel(' Amplitude Freq. Res.(dB)');
xlabel(' Digital Freq (rad)');
参考资料:数字信号处理-理论、算法与实现