Matlab平滑滤波&IIR音频滤波

鹅毛在路上了
优质创作者: 嵌入式与硬件开发技术领域
2023-05-03 16:02:13

目录

  • Matlab滑动平均滤波器设计
  • Matlab IIR滤波器设计
  • IIR音频滤波应用

Matlab滑动平均滤波器设计

clc,clear,close all;
%matlab滑动平均滤波器基础例程:
x = linspace(-pi,pi,100);
y = sin(2*pi*x);
n = randn(size(x));
t = y+n;
y1 = movmean(t,8);
figure
subplot(211);
plot(x,t,'LineWidth',1.1);
grid on
title('带噪声输入信号');
subplot(212);
plot(x,y1,'LineWidth',1.1);
grid on
title('平滑滤波后的信号');

% 滑动平均滤波器的幅频响应和相频响应
% 定义滤波器的阶数
Ns = [3,11,31];
for N = Ns
    % 定义采样频率和截止频率
    fs = 1000; % 采样频率为1000Hz
    fc = 50; % 截止频率为50Hz
    % 计算滤波器系数
    b = ones(1,N)/N; % 滑动平均滤波器的系数为1/N
    % 计算频率响应
    [H,f] = freqz(b,1,[],fs);
    % 计算幅频响应和相频响应
    mag = 20*log10(abs(H)); % 幅频响应
    phase = unwrap(angle(H))*(180/pi); % 相频响应,转化为角度制
    % 绘制幅频响应曲线
    figure;
    subplot(211)
    plot(f,mag,'LineWidth',1.1);
    grid on
    xlabel('频率 (Hz)');
    ylabel('幅值 (dB)');
    title(sprintf('N=%d的滑动平均滤波器幅频响应', N));
    % 绘制相频响应曲线
    subplot(212)
    plot(f,phase,'LineWidth',1.1);
    grid on
    xlabel('频率 (Hz)');
    ylabel('相位 (°)');
    title(sprintf('N=%d的滑动平均滤波器相频响应', N));
end

%读取"Don_Giovanni.wav"并进行滑动平均滤波测试
[music0,Fs] = audioread('Don_Giovanni.wav');
player = audioplayer(music0,Fs);
play(player)
% stop(player)  %停止播放
pause(abs(length(music0)/Fs))
disp('原始音频播放完毕')

music3 = movmean(music0,31);
player = audioplayer(music3,Fs); %N=3的效果
play(player)
pause(abs(length(music3)/Fs))
disp('N=3滑动平均滤波音频播放完毕')

music11 = movmean(music0,11);
player = audioplayer(music11,Fs); %N=11的效果
play(player)
pause(abs(length(music11)/Fs))
disp('N=11滑动平均滤波音频播放完毕')

music31 = movmean(music0,31);
player = audioplayer(music31,Fs); %N=31的效果
play(player)
pause(abs(length(music31)/Fs))
disp('N=31滑动平均滤波音频播放完毕')

img

img

img

Matlab IIR滤波器设计

clc,clear,close all;  
%**************************************************************
% Filtering a signal by a filter with a transfer function
%        B(z)    b0+b1z^{-1}+b2z^{-2}+ ... +bMz^{-M}
% H(z)= ----- = -------------------------------------
%        A(z)    a0+a1z^{-1}+a2z^{-2}+ ... +aNz^{-N}
%**************************************************************       
% General parameters 
Fe=44100;  %采样频率        
Te=1/Fe;   %周期        
N=Fe;      %序列长度             
t=0:Te:(N-1)*Te;        
f=-Fe/2:Fe/N:Fe/2-Fe/N; 
% 1st order IIR filter
a=-0.5;   
b0 =  1+a;
b1 =  0;
a0 =  1;
a1 =  a;
B1=[b0 b1];
A1=[a0 a1];
% Filtering a step signal
echelon=ones(1,N);
rep_echelon=filter(B1,A1,echelon);                                    
% Synthesis of a square signal
fc=250;                       
carre=1+square(2*pi*fc*t);    
rep_carre=filter(B1,A1,carre);

Npts=floor(Fe/fc);   
figure(1) 
%阶跃信号滤波
stem(t(1:Npts),echelon(1:Npts),'k');             
grid on; hold on;             
stem(t(1:Npts),rep_echelon(1:Npts),'r');                 
xlabel('Time (s)'); 
ylabel('Amplitude');
legend('Input signal','Output signal')
title('1st order filter : step signal filtering')
 
%方波信号滤波
figure(2)
stem(t(1:Npts),carre(1:Npts),'k');             
grid on; hold on;             
stem(t(1:Npts),rep_carre(1:Npts),'r');                 
xlabel('Time (s)'); 
ylabel('Amplitude');
legend('Input signal','Output signal')
title('1st order filter : square signal filtering')

img

img


% Filtering a signal by a filter with a transfer function
% B(z) b0+b1z^{-1}+b2z^{-2}+ ... +bMz^{-M}
% H(z)= ----- = -------------------------------------
% A(z) a0+a1z^{-1}+a2z^{-2}+ ... +aNz^{-N}
%**************************************************************
% 1st order IIR filter 一阶IIR滤波器

分子分母系数:
b0 = 0.5
b1 = 0
a0 = 1
a1 = -0.5

系统函数:
H(z) = (1/2+0)/(1-1/2z^-1)
= 1/2/(1-1/2
z^-1)
= 1/(2-z^-1)
= z/(2*z-1) = Y(z)/X(z) (传递函数:2y(n) - y(n-1) = x(n))

零极点分析
一阶零点:z = 0
一阶极点:z = 0.5

稳定性:系统稳定时的收敛域为z>0.5,收敛域包含单位圆,故该系统是稳定的。

IIR滤波器可以对信号进行低通、高通、带通或带阻滤波。
低通滤波器可以通过去除高频信号来平滑信号;
高通滤波器可以通过去除低频信号来突出高频信号;
带通和带阻滤波器可以选择性地保留或去除一定范围内的频率分量;
在本代码中,通过part2的波形图可以看出信号跳变部分变的更加平滑(变化更加缓慢),
使用的是一阶低通IIR滤波器(w=0(z=1)时,H(z)不为0.),其目的是平滑输入信号。
因此,滤波后的输出信号可能会比输入信号更平滑。

IIR音频滤波应用

clc,clear,close all;  
Fe=44100;  %采样频率        
Te=1/Fe;   %周期        
N=Fe;      %序列长度             
t=0:Te:(N-1)*Te;        
f=-Fe/2:Fe/N:Fe/2-Fe/N; 
% 2st order IIR filter 二阶IIR
a_1=-1.2;
a_2=.3;
b0 = 1+a_1+a_2;
b1 = 0;
b2 = 0;
a0 = 1;
a1 = a_1;
a2 = a_2;
B2=[b0 b1 b2];
A2=[a0 a1 a2];

% Filtering a step signal
echelon=ones(1,N);
rep_echelon=filter(B2,A2,echelon);                                    
% Synthesis of a square signal
fc=250;                       
carre=1+square(2*pi*fc*t);    
rep_carre=filter(B2,A2,carre);
% Sinusoidal signal
n = 1:N;
sin_signal=sin(2*pi*1000.*n./Fe);
rep_sin_signal=filter(B2,A2,sin_signal);

Npts=floor(Fe/fc);   
figure(1) 
%阶跃信号滤波
stem(t(1:Npts),echelon(1:Npts),'k');             
grid on; hold on;             
stem(t(1:Npts),rep_echelon(1:Npts),'r');                 
xlabel('Time (s)'); 
ylabel('Amplitude');
legend('Input signal','Output signal')
title('2st order filter : step signal filtering')
 
%方波信号滤波
figure(2)
stem(t(1:Npts),carre(1:Npts),'k');             
grid on; hold on;             
stem(t(1:Npts),rep_carre(1:Npts),'r');                 
xlabel('Time (s)'); 
ylabel('Amplitude');
legend('Input signal','Output signal')
title('2st order filter : square signal filtering') 

%正弦信号滤波
figure(3)
stem(t(1:Npts),sin_signal(1:Npts),'k');             
grid on; hold on; 
stem(t(1:Npts),rep_sin_signal(1:Npts),'r'); 
xlabel('Time (s)'); 
ylabel('Amplitude');
legend('Input signal','Output signal')
title('2st order filter : Sinusoidal signal filtering')

% 音乐信号滤波
[music0,Fs] = audioread('Don_Giovanni.wav');
player = audioplayer(music0,Fs);
play(player)
% stop(player)  %停止播放
pause(abs(length(music0)/Fs))
disp('原始音频播放完毕')

music1 = filter(B2,A2,music0);
player = audioplayer(music1,Fs);
play(player)
pause(abs(length(music1)/Fs))
disp('IIR滤波音频播放完毕')

figure(4)
subplot(2,1,1)
plot(music0)
grid on
title('原始音乐信号')
subplot(2,1,2)
plot(music1)
grid on
title('IIR滤波后的音乐信号')

img

img

...全文
316 回复 打赏 收藏 转发到动态 举报
写回复
用AI写文章
回复
切换为时间正序
请发表友善的回复…
发表回复

186

社区成员

发帖
与我相关
我的任务
社区描述
欢迎加入Matlab编程社区,支持各专业领域有关Matlab的博文、编程知识、经验分享,感谢每一份支持,您的鹅毛在路上了~
matlab 个人社区 山东省·潍坊市
社区管理员
  • 鹅毛在路上了
  • 数学建模加油站
  • MATLAB码农
加入社区
  • 近7日
  • 近30日
  • 至今
社区公告
暂无公告

试试用AI创作助手写篇文章吧