基于matlab通过多个FIR，IIR形成对violin信号进行均衡

基于《数字信号处理》这门课程的课程设计，特意写下这篇博客，仅供有需要的看官参考。

理论知识：

FIR和IIR滤波器是数字信号处理中常用的滤波器类型。设计带通滤波器时，需要考虑滤波器的通带、阻带、过渡带等参数，以及滤波器的阶数、截止频率等设计参数。

FIR滤波器是一种无限脉冲响应滤波器，其基本理论是通过对输入信号的加权和延迟来实现滤波。FIR滤波器的特点是稳定性好、易于设计和实现。

IIR滤波器是一种有限脉冲响应滤波器，其基本理论是通过对输入信号的加权和延迟以及对输出信号的反馈来实现滤波。IIR滤波器的特点是具有较高的性能和效率，但设计和实现相对复杂。

任务要求：设计三种不同增益的均衡器，分别增强低音、中音和高音

首先，我们需要通过matlab获取violin信号时域图和幅频特性：

clc;
clear;
[x0,fs]=audioread('violin.wav');%获取violin信号的波形数据以及采样率
N=length(x0);%获取音频信号长度
t=0:1/fs:(N-1)/fs;%得到时间轴
faxis=fs/N:fs/N:fs/2-fs/N;
%获取波形的时域图以及频谱图
figure
ffx0=fftshift(abs(fft(x0)));ffx=ffx0(N/2+2:N);
subplot(2,1,1);plot(t,x0,'k-','LineWidth',1.5);box off;
xlim([0 max(t)]);%ylim([-3 3]);
xlabel('Time (s)');ylabel('Amplitude');
subplot(2,1,2);plot(faxis,ffx,'k-','LineWidth',1.5);box off;

接着设计多个FIR滤波器实现均衡器：

clc;
clear;
[x0,fs]=audioread('violin.wav');
N=length(x0);
t=0:1/fs:(N-1)/fs;
faxis=fs/N:fs/N:fs/2-fs/N;
 
figure
ffx0=fftshift(abs(fft(x0)));ffx=ffx0(N/2+2:N);
subplot(2,1,1);plot(t,x0,'k-','LineWidth',1.5);box off;
xlim([0 max(t)]);%ylim([-3 3]);
xlabel('Time（s）');ylabel('Amplitude');
subplot(2,1,2);plot(faxis,ffx,'k-','LineWidth',1.5);box off;
 
%%  滤波器指标
rp=0.1;%rp表示通带最大衰减量
rs=0.001;%rs表示阻带最小衰减量
 
Ap=-20*log10(1-rp);
As=-20*log10(rs);%将其转换为dB
 %该频率是根据上述幅频特性曲线进行选取
wn1=[40 60 700 720];%设置具有过渡带宽为20的通带和阻带，[60,700]在该频率通过，
                    %小于40，大于720的不允许通过，下述9个与之相同
wn2=[740 760 1100 1120]; 
wn3=[1180 1200 1400 1420];
wn4=[1500 1520 1800 1820];  
wn5=[1840 1860 2000 2020]; 
wn6=[2040 2060 2340 2360];  
wn7=[2380 2400 2700 2720]; 
wn8=[2780 2800 3100 3120];  
wn9=[3270 3290 3600 3620];  
%% FIR 均衡器设计 用firpm
%%% 计算滤波器参数 firpmord
mags=[0 1 0];%描述带通
devs=[rs rp rs];%用于计算滤波器参数的值
[n1,f01,ao1,w1]=firpmord(wn1,mags,devs,fs);%计算滤波器参数
[n2,f02,ao2,w2]=firpmord(wn2,mags,devs,fs);
[n3,f03,ao3,w3]=firpmord(wn3,mags,devs,fs);
[n4,f04,ao4,w4]=firpmord(wn4,mags,devs,fs);
[n5,f05,ao5,w5]=firpmord(wn5,mags,devs,fs);
[n6,f06,ao6,w6]=firpmord(wn6,mags,devs,fs);
[n7,f07,ao7,w7]=firpmord(wn7,mags,devs,fs);
[n8,f08,ao8,w8]=firpmord(wn8,mags,devs,fs);
[n9,f09,ao9,w9]=firpmord(wn9,mags,devs,fs);
%%% 计算滤波器系数
hh1=firpm(n1,f01,ao1,w1);%计算滤波器系数
hh2=firpm(n2,f02,ao2,w2);
hh3=firpm(n3,f03,ao3,w3);
hh4=firpm(n4,f04,ao4,w4);
hh5=firpm(n5,f05,ao5,w5);
hh6=firpm(n6,f06,ao6,w6);
hh7=firpm(n7,f07,ao7,w7);
hh8=firpm(n8,f08,ao8,w8);
hh9=firpm(n9,f09,ao9,w9);
%%% 画出滤波器幅频特性曲线
[h1 W1]=freqz(hh1,1,1024);
[h2,W2]=freqz(hh2,1,1024);
[h3,W3]=freqz(hh3,1,1024);
[h4,W4]=freqz(hh4,1,1024);
[h5,W5]=freqz(hh5,1,1024);
[h6,W6]=freqz(hh6,1,1024);
[h7,W7]=freqz(hh7,1,1024);
[h8,W8]=freqz(hh8,1,1024);
[h9,W9]=freqz(hh9,1,1024);
figure; %绘图，均衡器幅频响应
loglog(W1,abs(h1),W2,abs(h2), W3,abs(h3),...
W4,abs(h4),W5,abs(h5),W6,abs(h6),W7,abs(h7),...
W8,abs(h8),W9,abs(h9));
xlabel('Frequency (Hz)r) ;ylabel(filter Gain');
axis([10 10^5   10^(-9) 1]);
%%%  对音乐信号进行均衡并画出频谱
y1=filter(hh1,1,x0);%设计滤波器对x0信号进行滤波
y2=filter(hh2,1,x0);
y3=filter(hh3,1,x0);
y4=filter(hh4,1,x0);
y5=filter(hh5,1,x0);
y6=filter(hh6,1,x0);
y7=filter(hh7,1,x0);
y8=filter(hh8,1,x0);
y9=filter(hh9,1,x0);
g1= 10;g2 = 10;g3 = 10;g4 = 0;
g5 = 0;g6 = 0;g7 = 0;g8 = 0;g9 = 0;
y=y1.*g1+y2.*g2+y3.*g3+y4.*g4+y5.*g5+y6.*g6+y7.*g7+y8.*g8+y9.*g9+x0;
%对前三个信号即低音进行均衡则选取g1,g2,g3设置为10，依次向后推理获得中音和高音均衡
%对均衡信号画出时域和频域波形图
figure
ffy=fftshift(abs(fft(y)));ffy0=ffy(N/2+2:N);
subplot(2,1,1);plot(t,y,'k-','LineWidth',1.5);box off;
xlim([0 max(t)]);%ylim([-3 3]);
xlabel('Time（s）');ylabel('Amplitude');
subplot(2,1,2);plot(faxis,ffy0,'k-','LineWidth',1.5);box off;
%%% 听效果并存储
sound(y);

该均衡器获得的幅频特性曲线如下图：

对低音进行均衡后信号的时域频域波形图如下：

 

对中音进行均衡后信号的时域频域波形图如下：

 

对高音进行均衡后信号的时域频域波形图如下：

基于matlab的多个IIR设计实现均衡器：

clc;
clear;
[x0,fs]=audioread('violin.wav');
N=length(x0);
t=0:1/fs:(N-1)/fs;
faxis=fs/N:fs/N:fs/2-fs/N;

figure
ffx0=fftshift(abs(fft(x0)));ffx=ffx0(N/2+2:N);
subplot(2,1,1);plot(t,x0,'k-','LineWidth',1.5);box off;
xlim([0 max(t)]);%ylim([-3 3]);
xlabel('Time (s)');ylabel('Amplitude');
subplot(2,1,2);plot(faxis,ffx,'k-','LineWidth',1.5);box off;

%% 滤波器指标
Ap=1;
As=20;


wp1=[60 700]; ws1=[40 720];
wp2=[760 1100]; ws2=[740 1120];
wp3=[1200 1400]; ws3=[1180 1420];
wp4=[1520 1800]; ws4=[1500 1820];
wp5=[1860 2000]; ws5=[1840 2020];
wp6=[2060 2340]; ws6=[2040 2360];
wp7=[2400 2700]; ws7=[2380 2720];
wp8=[2800 3100]; ws8=[2780 3120];
wp9=[3290 3600]; ws9=[3270 3620];
wp10=[3660 6000];ws10=[3640 6020];

%% IIR 均衡器设计 ellip
wpI1=2*wp1/fs;wsI1=2*ws1/fs;%将频率进行归一化
wpI2=2*wp2/fs;wsI2=2*ws2/fs;
wpI3=2*wp3/fs;wsI3=2*ws3/fs;
wpI4=2*wp4/fs;wsI4=2*ws4/fs;
wpI5=2*wp5/fs;wsI5=2*ws5/fs;
wpI6=2*wp6/fs;wsI6=2*ws6/fs;
wpI7=2*wp7/fs;wsI7=2*ws7/fs;
wpI8=2*wp8/fs;wsI8=2*ws8/fs;
wpI9=2*wp9/fs;wsI9=2*ws9/fs;
wpI10=2*wp10/fs;wsI10=2*ws10/fs;

%%% 计算滤波器参数 ellipord
[n1,wn1]=ellipord(wpI1,wsI1,Ap,As);Niir1=n1;
[n2,wn2]=ellipord(wpI2,wsI2,Ap,As);Niir2=n2;
[n3,wn3]=ellipord(wpI3,wsI3,Ap,As);Niir3=n3;
[n4,wn4]=ellipord(wpI4,wsI4,Ap,As);Niir4=n4;
[n5,wn5]=ellipord(wpI5,wsI5,Ap,As);Niir5=n5;
[n6,wn6]=ellipord(wpI6,wsI6,Ap,As);Niir6=n6;
[n7,wn7]=ellipord(wpI7,wsI7,Ap,As);Niir7=n7;
[n8,wn8]=ellipord(wpI8,wsI8,Ap,As);Niir8=n8;
[n9,wn9]=ellipord(wpI9,wsI9,Ap,As);Niir9=n9;
[n10,wn10]=ellipord(wpI10,wsI10,Ap,As);Niir10=n10;

%%% 计算滤波器系数
[b1,a1]=ellip(n1,Ap,As,wn1);
[b2,a2]=ellip(n2,Ap,As,wn2);
[b3,a3]=ellip(n3,Ap,As,wn3);
[b4,a4]=ellip(n4,Ap,As,wn4);
[b5,a5]=ellip(n5,Ap,As,wn5);
[b6,a6]=ellip(n6,Ap,As,wn6);
[b7,a7]=ellip(n7,Ap,As,wn7);
[b8,a8]=ellip(n8,Ap,As,wn8);
[b9,a9]=ellip(n9,Ap,As,wn9);
[b10,a10]=ellip(n10,Ap,As,wn10);

%%% 画出滤波器幅频特性曲线
[h1,W1]=freqz(b1,a1);
[h2,W2]=freqz(b2,a2);
[h3,W3]=freqz(b3,a3);
[h4,W4]=freqz(b4,a4);
[h5,W5]=freqz(b5,a5);
[h6,W6]=freqz(b6,a6);
[h7,W7]=freqz(b7,a7);
[h8,W8]=freqz(b8,a8);
[h9,W9]=freqz(b9,a9);
[h10,W10]=freqz(b10,a10);

figure; %绘图，均衡器幅频响应
loglog(W1,abs(h1),W2,abs(h2), W3,abs(h3),...
W4,abs(h4),W5,abs(h5),W6,abs(h6),W7,abs(h7),...
W8,abs(h8),W9,abs(h9),W10,abs(h10));
xlabel('Frequency (Hz)'); ylabel('Filter Gain');
axis([10 10^5 10^-9 1]);

%%% 对音乐信号进行均衡并画出频谱
y1=filter(b1,a1,x0);
y2=filter(b2,a2,x0);
y3=filter(b3,a3,x0);
y4=filter(b4,a4,x0);
y5=filter(b5,a5,x0);
y6=filter(b6,a6,x0);
y7=filter(b7,a7,x0);
y8=filter(b8,a8,x0);
y9=filter(b9,a9,x0);
y10=filter(b10,a10,x0);
g1=0; g2=0; g3=0; g4=0;
g5=10; g6=10; g7=10; g8=0; g9=0; g10=0;
y = y1.*g1 + y2.*g2 + y3.*g3 + y4.*g4 + y5.*g5 +...
y6.*g6 + y7.*g7 + y8.*g8 + y9.*g9 + y10.*g10 + x0;

figure
ffy=fftshift(abs(fft(y)));ffy0=ffy(N/2+2:N);
subplot(2,1,1);plot(t,y,'k-','LineWidth',1.5);box off;
xlim([0 max(t)]);%ylim([-3 3]);
xlabel('Time（s）');ylabel('Amplitude');
subplot(2,1,2);plot(faxis,ffy0,'k-','LineWidth',1.5);box off;
%%% 听效果并存储
sound(y,fs);

该均衡器获得的幅频特性曲线如下图：

IIR设计均衡器对低音进行增强后均衡处理后信号的时域，频域图如下：

IIR设计均衡器对中音进行增强后均衡处理后信号的时域，频域图如下：

IIR设计均衡器对高音进行增强后均衡处理后信号的时域，频域图如下：

上述代码有点呆（很早之前写的代码），各位优化代码可通过数组或for循环进行简化，

另外，若存在错误，请各位多多指正！