matlab设计FIR滤波器

方法1:通过fir1()函数进行设计

        B = fir1(N,Wn)设计FIR低通滤波器，返回的滤波器参数保存在长度为N+1的数组B中。Wn为归一化截止频率，范围为0~1。截止频率用于区分过渡带和阻带。1处对应的是采样频率的一半。滤波器系数B是实的且有线性相位。

B=fir1(N,Wn,'high'); %表示N阶高通滤波器   B=fir1(N,Wn,'low'); %表示N阶低通滤波器;

如果Wn=[W1, W2]; B=fir1(N,Wn,'bandpass');  %N阶带通滤波器 B=fir1(N,Wn,'stop'); %N阶带阻滤波器。

freqz(B)；%绘制滤波器的频域波形

sin_signal_filter = filter(filter_lowpass,1,sin_signal); %滤波函数，参数1为滤波器系数，参数2为行数（这里没有使用矩阵）故为1，参数3为输入信号。

测试代码1：

%%%%%%%%%%%%%%%%%%% 低通滤波器仿真程序 %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% File: fir_sim.m   %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%% Date:2023-5-5 Author:侯永明 %%%%%%%%%%%%%
%%%%%*******************程序说明**********************%%%%%%%
%本程序完成了对于单音信号和双音信号的低通滤波仿真
%版本信息：MATLAB 2017a

%%%%%*******************程序主体**********************%%%%%%%
clear clc
close all;

%%%%参数设置
time=1:1000; %信号长度
fre1 =50;    %单音信号1的频率
fre2 =500;   %单音信号2的频率
fs=8000;     %采样频率

sin_signal1 = sin(2*pi*fre1*time/fs);
sin_signal2 = sin(2*pi*fre2*time/fs);
sin_signal = sin_signal1+sin_signal2;

%%%%%观察单音信号和双音信号的时域图和频域图
figure(1)
plot(sin_signal1);
title('单音信号1的时域图');

figure(2)
plot(sin_signal2);
title('单音信号2的时域图');

figure(3)
plot(sin_signal);
title('双音信号的时域图');

figure(4)
plot(abs(fft(sin_signal1)));
title('单音信号1的频域图');

figure(5)
plot(abs(fft(sin_signal2)));
title('单音信号2的频域图');

figure(6)
plot(abs(fft(sin_signal)));
title('双音信号的频域图');

%%%%%设计低通滤波器
filter_lowpass = fir1(128,550/8000); %归一化截止频率的公式：w=fc*2/fs
                                     %128为滤波阶数 fir1是滤波器的构造函数
figure(7)
freqz(filter_lowpass);               %freqz是求频域图形的函数
title('滤波器频域响应图');

%%%%%对双音信号进行滤波
sin_signal_filter = filter(filter_lowpass,1,sin_signal);

%%%%%观察滤波之后的时域图和频域图
figure(8)
plot(sin_signal_filter);
title('滤波之后信号的时域图');

figure(9)
plot(abs(fft(sin_signal_filter)));
title('滤波之后信号的频域图');         
%%%%%*******************程序结束**********************%%%%%%%

测试代码2:

clear all;
close all;
clc;

%滤波器长度
N=41;
%采样频率
fs=2000;

%各种滤波器的频率特征
fc_lpf = 200;
fc_hpf = 200;
fp_bandpass = [200 400];
fc_stop = [200 400];

%以采样频率的一半，对频率进行归一化
wn_lpf = fc_lpf*2/fs;
wn_hpf = fc_hpf*2/fs;
wn_bandpass = fp_bandpass*2/fs;
wn_stop = fc_stop*2/fs;

%采用fir1函数设计fir滤波器
b_lpf = fir1(N-1,wn_lpf);
b_hpf = fir1(N-1,wn_hpf,'high');
b_bandpass = fir1(N-1,wn_bandpass, 'bandpass');
b_stop = fir1(N-1,wn_stop,'stop');

%求幅频响应
m_lpf = 20*log(abs(fft(b_lpf)))/log(10);
m_hpf = 20*log(abs(fft(b_hpf)))/log(10);
m_bandpass = 20*log(abs(fft(b_bandpass)))/log(10);
m_stop = 20*log(abs(fft(b_stop)))/log(10);

%设置频率响应单位的横坐标为Hz
x_f = 0:(fs/length(m_lpf)):fs/2;        %???

%绘制单位脉冲响应
subplot(4,2,1); stem(b_lpf);xlabel('n');ylabel('h(n)');legend('lpf'); %stem绘制b_lpf的离散序列图 legend函数为数据创建图例
subplot(4,2,3); stem(b_hpf);xlabel('n');ylabel('h(n)');legend('hpf');
subplot(4,2,5); stem(b_bandpass);xlabel('n');ylabel('h(n)');legend('bandpass');
subplot(4,2,7); stem(b_stop);xlabel('n');ylabel('h(n)');legend('stop');

%绘制幅频响应
subplot(4,2,2); plot(x_f,m_lpf(1:length(x_f))); xlabel('频率(Hz)');ylabel('幅度(dB)','fontsize',8); legend('lpf');
subplot(4,2,4); plot(x_f,m_hpf(1:length(x_f))); xlabel('频率(Hz)');ylabel('幅度(dB)','fontsize',8); legend('hpf');
subplot(4,2,6); plot(x_f,m_bandpass(1:length(x_f))); xlabel('频率(Hz)');ylabel('幅度(dB)','fontsize',8); legend('bandpass');
subplot(4,2,8); plot(x_f,m_stop(1:length(x_f))); xlabel('频率(Hz)');ylabel('幅度(dB)','fontsize',8); legend('stop');

测试代码3：

%%低通，过渡带为1~1.5KHz,采样频率8kHz,通带最大波纹0.01，阻带最大波纹0.05

%获取凯赛窗参数，后续fir1要设计凯赛窗的低通滤波器

fs=8000;           %采样频率
fc=[1000,1500];    %过渡带; 1000-1500过渡带. 通带:0~（1000*2/fs）,阻带:(1500*2/fs)
mag=[1 0];         %理想幅值 0~（1000*2/fs）是1，1500*2/fs~1是0
dev=[0.01,0.05];   %波纹
[n,wn,beta,ftype] = kaiserord(fc,mag,dev,fs); %获取凯赛窗参数，根据据过渡带幅值返回ftype类型

%fir1设计凯赛窗和汉明窗滤波器
h_kaiser=fir1(n,wn,ftype,kaiser(n+1,beta));
h_hamm=fir1(n,fc(2)*2/fs);

%设计最优滤波器
fpm=[0 fc(1)*2/fs fc(2)*2/fs 1]; %firmp频段向量，归一化[0 1]
magpm =[1 1 0 0]; %firpm幅值向量
h_pm= firpm(n,fpm,magpm); %设计低通滤波器

%幅值
m_kaiser=20*log(abs(fft(h_kaiser,1024)))/log(10);
m_hamm=20*log(abs(fft(h_hamm,1024)))/log(10);
m_pm=20*log(abs(fft(h_pm,1024)));

%设置幅频响应的横坐标为Hz
x_f=[0:1:length(m_kaiser)/2]*fs/length(m_kaiser);

%只显示正频率部分
m1=m_kaiser(1:length(x_f));
m2=m_hamm(1:length(x_f));
m3=m_pm(1:length(x_f));

%绘制幅频曲线
plot(x_f,m1, '-',x_f,m2,'-.',x_f,m3,'--');
xlabel('频率(Hz)');
ylabel('幅度(dB)');
legend('凯赛窗','汉明窗','最优滤波');
grid;

  2.方法2：通过fdatool（filter designer）来实现滤波器的设计。

在工具箱弹窗中直接进行参数的选择和设置，可以选择滤波器类型、设计方法、阶数、通带阻带频率、采样频率及幅值。

 设置完成之后，点击下图的文件->生成matlab代码->滤波器设计函数进行.m文件保存

 可以生成 .m文件，需要注意把第一行函数头手动注释掉，代码如下：

%function Hd = ditong2
%DITONG2 Returns a discrete-time filter object.

% MATLAB Code
% Generated by MATLAB(R) 9.2 and the Signal Processing Toolbox 7.4.
% Generated on: 05-May-2023 14:26:21

% Equiripple Lowpass filter designed using the FIRPM function.

% All frequency values are in Hz.
Fs = 1000;  % Sampling Frequency

Fpass = 150;             % Passband Frequency                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            
Fstop = 200;             % Stopband Frequency
Dpass = 0.057501127785;  % Passband Ripple
Dstop = 0.0001;          % Stopband Attenuation
dens  = 20;              % Density Factor

% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fpass, Fstop]/(Fs/2), [1 0], [Dpass, Dstop]);

% Calculate the coefficients using the FIRPM function.
b  = firpm(N, Fo, Ao, W, {dens});
Hd = dfilt.dffir(b);

% [EOF]

 Run一下生成的滤波器函数.m文件，会生成滤波器系数矩阵Hd，然后就可以调用filter函数对输入信号进行滤波了。

%%
y1 = rand(50,1);
subplot(4,1,1);
plot(y1);
title("滤波前波形");
subplot(4,1,2);
plot(abs(fft(y1)));
title("滤波前频谱");
%%
y2 = filter(Hd,y1);
subplot(4,1,3);
plot(abs(fft(y2)));
title("滤波后频谱");
subplot(4,1,4);
plot(y2);
title("滤波后波形");