7-band EQ均衡器算法

函数

fs=44.1*10^3;
% fs=48*10^3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EQ_band1
gain1=-10 ;bw1=3000;fc1=10000;
k1=10^(gain1/20);acj1=(1-sin(2*pi*bw1/fs))/cos(2*pi*bw1/fs);
bcj1=cos(2*pi*fc1/fs);
% b1(1)=((1+k1)+(1-k1)*acj1)/2;
% b1(2)=-(1+acj1)*bcj1;
% b1(3)=((1-k1)+(1+k1)*acj1)/2;
% a1(1)=1;
% a1(2)=-(1+acj1)*bcj1;
% a1(3)=acj1;
b1(1)=((1+k1)+(1-k1)*acj1)/2;
b1(2)=-(1+acj1)*bcj1  /2;
b1(3)=((1-k1)+(1+k1)*acj1)/2;
a1(1)=1;
a1(2)=-(1+acj1)*bcj1   /(-2);
a1(3)=acj1   /(-1);
%b1(1)=1;b1(2)=0;b1(3)=0;
%a1(1)=1;a1(2)=0;a1(3)=0;
%[ampo1,freqo1]=freqz(b1,a1);
%plot(freqo1*fs*0.5/pi,20*log10(abs(ampo1)));
%set(gca,'Xscale','log');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EQ_band2
gain2=7;bw2=200;fc2=500;
k2=10^(gain2/20);acj2=(1-sin(2*pi*bw2/fs))/cos(2*pi*bw2/fs);
bcj2=cos(2*pi*fc2/fs);
b2(1)=((1+k2)+(1-k2)*acj2)/2;
b2(2)=-(1+acj2)*bcj2;
b2(3)=((1-k2)+(1+k2)*acj2)/2;
a2(1)=1;
a2(2)=-(1+acj2)*bcj2;
a2(3)=acj2;
%b2(1)=1;b2(2)=0;b2(3)=0;
%a2(1)=1;a2(2)=0;a2(3)=0;
%[ampo2,freqo2]=freqz(b2,a2);
%plot(freqo2*fs*0.5/pi,20*log10(abs(ampo2)));
%set(gca,'Xscale','log');
[b_all1,a_all1]=series(b1,a1,b2,a2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EQ_band3
gain3=8;bw3=250;fc3=1000;
k3=10^(gain3/20);acj3=(1-sin(2*pi*bw3/fs))/cos(2*pi*bw3/fs);
bcj3=cos(2*pi*fc3/fs);
b3(1)=((1+k3)+(1-k3)*acj3)/2;
b3(2)=-(1+acj3)*bcj3;
b3(3)=((1-k3)+(1+k3)*acj3)/2;
a3(1)=1;
a3(2)=-(1+acj3)*bcj3;
a3(3)=acj3;
%b3(1)=1;b3(2)=0;b3(3)=0;
%a3(1)=1;a3(2)=0;a3(3)=0;
%[ampo3,freqo3]=freqz(b3,a3);
%plot(freqo3*fs*0.5/pi,20*log10(abs(ampo3)));
%set(gca,'Xscale','log');
[b_all2,a_all2]=series(b_all1,a_all1,b3,a3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EQ_band4
gain4=9;bw4=300;fc4=2000;
k4=10^(gain4/20);acj4=(1-sin(2*pi*bw4/fs))/cos(2*pi*bw4/fs);
bcj4=cos(2*pi*fc4/fs);
b4(1)=((1+k4)+(1-k4)*acj4)/2;
b4(2)=-(1+acj4)*bcj4;
b4(3)=((1-k4)+(1+k4)*acj4)/2;
a4(1)=1;
a4(2)=-(1+acj4)*bcj4;
a4(3)=acj4;
%b4(1)=1;b4(2)=0;b4(3)=0;
%a4(1)=1;a4(2)=0;a4(3)=0;
%[ampo4,freqo4]=freqz(b4,a4);
%plot(freqo4*fs*0.5/pi,20*log10(abs(ampo4)));
%set(gca,'Xscale','log');
[b_all3,a_all3]=series(b_all2,a_all2,b4,a4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EQ_band5
gain5=10;bw5=350;fc5=3000;
k5=10^(gain5/20);acj5=(1-sin(2*pi*bw5/fs))/cos(2*pi*bw5/fs);
bcj5=cos(2*pi*fc5/fs);
b5(1)=((1+k5)+(1-k5)*acj5)/2;
b5(2)=-(1+acj5)*bcj5;
b5(3)=((1-k5)+(1+k5)*acj5)/2;
a5(1)=1;
a5(2)=-(1+acj5)*bcj5;
a5(3)=acj5;
%b5(1)=1;b5(2)=0;b5(3)=0;
%a5(1)=1;a5(2)=0;a5(3)=0;
%[ampo5,freqo5]=freqz(b5,a5);
%plot(freqo5*fs*0.5/pi,20*log10(abs(ampo5)));
%set(gca,'Xscale','log');
[b_all4,a_all4]=series(b_all3,a_all3,b5,a5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EQ_band6 lpf
gain6=11;fc6=8000;
k6=10^(gain6/20);acj6=(1-sin(2*pi*fc6/fs))/cos(2*pi*fc6/fs);
b6(1)=((1+k6)+(1-k6)*acj6)/2;
b6(2)=((-1+k6)-(1+k6)*acj6)/2;
b6(3)=0 ;
a6(1)=1;
a6(2)=-acj6;
a6(3)=0 ;
%b6(1)=1;b6(2)=0;b6(3)=0;
%a6(1)=1;a6(2)=0;a6(3)=0;
%[ampo6,freqo6]=freqz(b6,a6);
%plot(freqo6*fs*0.5/pi,20*log10(abs(ampo6)));
%set(gca,'Xscale','log');
[b_all5,a_all5]=series(b_all4,a_all4,b6,a6);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EQ_band7 hpf
gain7=12;fc7=60;
k7=10^(gain7/20);acj7=(1-sin(2*pi*fc7/fs))/cos(2*pi*fc7/fs);
b7(1)=((1+k7)-(1-k7)*acj7)/2;
b7(2)=((1-k7)-(1+k7)*acj7)/2;
b7(3)=0 ;
a7(1)=1;
a7(2)=-acj7;
a7(3)=0 ;
%b7(1)=1;b7(2)=0;b7(3)=0;
%a7(1)=1;a7(2)=0;a7(3)=0;
%[ampo7,freqo7]=freqz(b7,a7);
%plot(freqo7*fs*0.5/pi,20*log10(abs(ampo7)));
%set(gca,'Xscale','log');
[b_all6,a_all6]=series(b_all5,a_all5,b7,a7);
%[ampo8,freqo8]=freqz(b_all6,a_all6);
%plot(freqo8*fs*0.5/pi,20*log10(abs(ampo8)));
%set(gca,'Xscale','log');

主程序

clear;clc;
[y_org,fs,bits]=wavread('in.wav');
eq_calc;
y=filter(b_all6,a_all6,y_org(:,1));
wavwrite(y,fs,'out.wav');