数字滤波器设计

这又是一篇懒人福利博客嘿嘿嘿😝

设计步骤概览

• 使用designfilt函数设计/修改一个滤波器

  df = designfilt(resp,Name,Value);   %设计滤波器
  designfilt(df); %修改已存在滤波器

• 使用filter函数对指定信号序列进行滤波测试

• 使用fvtool工具箱可视化设计的滤波器

designfilt函数使用方法

df = designfilt(resp,Name,Value);

第一个参数为滤波器类型，然后是Name+Value的属性配置键值对（可以有很多组）

主要配置的重要属性有：滤波器类型、频率约束、幅度约束、滤波器阶数、设计算法、采样率

• 滤波器类型resp

  • lowpassfir
  • lowpassiir
  • highpassfir
  • highpassiir
  • bandpassfir
  • bandpassiir
  • bandstopfir
  • bandstopiir
  • differentiatorfir
  • hilbertfir
  • arbmagfir

• 属性配置键值对Name+Value

  在编写代码的时候键值对都是用逗号隔开噢

  属性	Name键	Value值类型
  滤波器阶数	FilterOrder	正整数
  分子阶数	NumeratorOrder	正整数
  分母阶数	DenominatorOrder	正整数
  通带频率	PassbandFrequency	正数
  阻带频率	StopbandFrequency	正数
  6dB截止频率	CutoffFrequency	正数
  3dB截止频率	HalfPowerFrequency	正数
  通带纹波	PassbandRipple	正数
  阻带衰减	StopbandAttenuation	正数
  采样率	SampleRate	正数

  • 设计方式DesignMethod
    • butter：具有通带内最平稳特性的Butterworth滤波器，但是衰减速度慢
    • cheby1：通带等纹波、阻带最平缓的切比雪夫Ⅰ型滤波器
    • cheby2：通带最平缓、阻带等纹波的切比雪夫Ⅱ型滤波器
    • cls：最小二乘约束滤波器设计FIR滤波器
    • ellip：通带阻带都等纹波的IIR椭圆滤波器
    • equiripple：使用Parks-McClellan算法设计FIR滤波器，可最小化各个带上的最大纹波
    • freqsamp：通过均匀采样频率响应并进行傅里叶反变换设计任意幅度响应的FIR滤波器
    • kaiserwin：使用Kaiser窗函数法
    • ls：最小二乘法设计FIR滤波器
    • maxflat：通带最平坦
    • window：先使用最小二乘法计算滤波器需要的参数然后使用你指定的窗函数平滑冲击响应曲线

设计实例

设计一个带阻滤波器滤掉指定几个频点

clear
%查看音频数据信息
sourceWaveInf = audioinfo('SunshineSquare.wav');
%读取原声音数据
[audio_data,fs] = audioread('SunshineSquare.wav');

%绘制音频波形
figure
subplot(2,2,1)
t = 0:seconds(1/fs):seconds(sourceWaveInf.Duration);    %时间轴数组
t = t(1:end-1);
plot(t,audio_data)  %画图
grid on
xlabel('Time')  %横轴时间
ylabel('Audio Signal')  %纵轴幅度
title('原音频信号时域波形')  %图表名

%频谱分析
if mod( sourceWaveInf.TotalSamples ,2 ) == 1
    %频谱需要对称的长度，所以要偶数
    sourceWavaLength = sourceWaveInf.TotalSamples - 1;
else
    sourceWavaLength = sourceWaveInf.TotalSamples;
end
subplot(2,2,2)
audioSpectrum = fft(audio_data);    %快速傅里叶变换计算
P2 = abs(audioSpectrum/sourceWavaLength);   %计算双侧频谱
P1 = P2(1:sourceWavaLength/2+1);    %取出单边频谱
P1(2 : end-1) = 2*P1(2 : end-1);    %双边合成单边，就是翻倍了
f = fs*(0:(sourceWavaLength/2))/sourceWavaLength;   %频率轴数组
plot(f,P1)  %画图
grid on
title('滤波前频谱')
xlabel('频率/Hz')
ylabel('FFT')

%设计滤波器
bsFilt_1 = designfilt('bandstopiir','FilterOrder',10,...
'HalfPowerFrequency1',1500,'HalfPowerFrequency2',1650,...
'SampleRate',fs);   %IIR滤波器，滤除1575Hz频率噪点
audioClear_1 = filter(bsFilt_1,audio_data); %滤波

bsFilt_2 = designfilt('bandstopiir','FilterOrder',10,...
'HalfPowerFrequency1',3100,'HalfPowerFrequency2',3200,...
'SampleRate',fs);   %IIR滤波器，滤除3150Hz频率噪点
audioClear_2 = filter(bsFilt_2,audioClear_1); %滤波

% bsFilt_3 = designfilt('bandstopiir','FilterOrder',10,...
% 'HalfPowerFrequency1',4700,'HalfPowerFrequency2',4750,...
% 'SampleRate',fs);   %IIR滤波器，滤除4725频率噪点
bsFilt_3 = designfilt('lowpassiir','FilterOrder',8, ...
        'PassbandFrequency',4000,'PassbandRipple',0.2, ...
        'SampleRate',11025);
audioClear = filter(bsFilt_3,audioClear_2); %滤波


%绘制滤波之后的音频信号
subplot(2,2,3)
plot(t,audioClear)
grid on
xlabel('Time')
ylabel('Audio Signal')
title('滤波后的信号')

%绘制滤波后的频谱
subplot(2,2,4)
audioSpectrum = fft(audioClear);    %快速傅里叶变换计算
P2 = abs(audioSpectrum/sourceWavaLength);   %计算双侧频谱
P1 = P2(1:sourceWavaLength/2+1);    %取出单边频谱
P1(2 : end-1) = 2*P1(2 : end-1);    %双边合成单边，就是翻倍了
f = fs*(0:(sourceWavaLength/2))/sourceWavaLength;   %频率轴数组
plot(f,P1)  %画图
title('滤波后频谱')
xlabel('频率/Hz')
ylabel('FFT')
grid on

%存储滤波后的音频数据
audiowrite('sunshineSquare_Clear.wav',audioClear,fs);
%播放滤波后的声音
sound(audioClear,fs);
%可视化设计的滤波器
% fvtool(bsFilt_1);

这里就不放图啦🎈笔者把这个小DEMO上传到GitHub仓库Matlab脚本库里面了欢迎大家白嫖

参考链接

• matlab官网-数字滤波器
• matlab官网-designfilt函数说明

如果觉得官网的全英文难以接受，欢迎联系笔者噢😘