%file:magnitude_response_equalizer.m;all-pass filter based magnitude response equalizer demo
clear,figure(1),figure(2),figure(3),figure(4),figure(5),figure(6),figure(7),clf
set(gcf,'doublebuffer','on')
w=0:0.002:2*pi;
UnitCircle=exp(j*w);
fs = 44100;
faxis=fs*w/(2*pi);
vBf=[100 350 1500 5500 9100];
vf0=[62.5 250 1000 4000 16000];
choicelist=str2mat('There are 5 channls in the frequency band 20Hz-20KHz:',...
' f01=62.5 Hz,',...
' f02=250 Hz,',...
' f03=1000 Hz,',...
' f04=4000 Hz,',...
' f05=16000 Hz,',...
' Please,make your choice from the list:',...
'1 - Press 1 for continue without any corrections.',...
'2 - Press 2 for continue with corrections.',...
'3 - Exit.');
maxchoice=size(choicelist,1)-1;
exitFlag=0;
while ~exitFlag
choice=0;
while isempty(choice)|(choice<1)|(choice>maxchoice)
disp(choicelist)
choice=input('Please,make your choice:');
end
switch choice
case 1
information=str2mat('All channls are without corrections.',...
'There is equality frequency response for the all frequency band.');
disp(information)
for n=1:5
K=1;
Bf=vBf(1,n);
Bw=2*pi*Bf/fs;
w0n=2*pi*vf0(1,n)/fs;
k1n=-cos(w0n);k2n=(1-tan(Bw/2))/(1+tan(Bw/2));
disp(['f0=' num2str(vf0(1,n)) ',k1=' num2str(k1n) ',k2=' num2str(k2n) ',K=' num2str(K)])
a=[1 k1n*(1+k2n) k2n];b=fliplr(a);
poles=roots(a);z=roots(b);
figure(n)
subplot(131),hold off,plot(UnitCircle,'g'),hold on,grid
plot(poles,'xr'),hz=plot(z,'o');set(hz,'MarkerSize',3)
xlabel('real(z)'),ylabel('imag(z)'),title('pole-zero plot')
H0=freqz(b,a,w);phi=unwrap(angle(H0));
subplot(132),semilogx(faxis,phi),grid,axis([0 fs/2 -2*pi 0])
xlabel('f,Hz'),ylabel('phase,rad'),title('phase response')
Fn=0.5*(1+K+(1-K)*H0);
subplot(133),semilogx(faxis,20*log10(abs(Fn)+eps))
grid,axis([0 fs/2 -10 10])
xlabel('f,Hz'),ylabel('K(f), dB'),title('frequency response')
if n==1
F1=Fn;k11=k1n;k21=k2n;
elseif n==2
F2=Fn;k12=k1n;k22=k2n;
elseif n==3
F3=Fn;k13=k1n;k23=k2n;
elseif n==4
F4=Fn;k14=k1n;k24=k2n;
else
F5=Fn;k15=k1n;k25=k2n;
end
end
F=F1.*F2.*F3.*F4.*F5;
figure(6)
semilogx(faxis,20*log10(abs(F)+eps))
grid,axis([0 fs/2 -10 10])
xlabel('f,Hz'),ylabel('K(f), dB'),title('total frequency response')
information1=str2mat('You can test the magnitide response equalizer',...
'with some real voice signal.',...
'Please,make your choice:',...
'1-Press 1 for testing the equalizer',...
'2-Press 2 for Exit.');
disp(information1)
choice1=input('Please,make your choice:');
switch choice1
case 1
load handel
sos=[0.5*((1+k21)+K*(1-k21)) k11*(1+k21) 0.5*((1+k21)-K*(1-k21)) 1 k11*(1+k21) k21
0.5*((1+k22)+K*(1-k22)) k12*(1+k22) 0.5*((1+k22)-K*(1-k22)) 1 k12*(1+k22) k22
0.5*((1+k23)+K*(1-k23)) k13*(1+k23) 0.5*((1+k23)-K*(1-k23)) 1 k13*(1+k23) k23
0.5*((1+k24)+K*(1-k24)) k14*(1+k24) 0.5*((1+k24)-K*(1-k24)) 1 k14*(1+k24) k24
0.5*((1+k25)+K*(1-k25)) k15*(1+k25) 0.5*((1+k25)-K*(1-k25)) 1 k15*(1+k25) k25];
figure(7)
subplot(121),psd(y);
title('Power Spectral Density Estimate via Wellch for the input signal')
y1=sosfilt(sos,y);
subplot(122),psd(y1);
title('Power Spectral Density Estimate via Wellch for the output signal')
information3=str2mat('You can hear the input and output signal.',...
'Please,make your choice:',...
'1-Press 1 for hearing the signals',...
'2-Press 2 for Exit.');
disp(information3)
choice2=input('Please,make your choice:');
switch choice2
case 1
wavplay(y)
wavplay(y1)
exitFlag=1;
case 2,exitFlag=1;
end
case 2,exitFlag=1;
end
case 2
information4=str2mat('All channls are able to be corrected',...
'with arbitrarily level.');
disp(information4)
s2='Please,enter the correction for ';
s3='channl f0=62.5 Hz';
s4='channl f0=250 Hz';
s5='channl f0=1000 Hz';
s6='channl f0=4000 Hz';
s7='channl f0=16000 Hz';
s8=' in dB:';
str2=[s2 s3 s8];
for n=1:5
dB=input(str2);
Bf=vBf(1,n);
Bw=2*pi*Bf/fs;
w0n=2*pi*vf0(1,n)/fs;
k1n=-cos(w0n);k2n=(1-tan(Bw/2))/(1+tan(Bw/2));
a=[1 k1n*(1+k2n) k2n];b=fliplr(a);
poles=roots(a);z=roots(b);
figure(n)
subplot(131),hold off,plot(UnitCircle,'g'),hold on,grid
plot(poles,'xr'),hz=plot(z,'o');set(hz,'MarkerSize',3)
xlabel('real(z)'),ylabel('imag(z)'),title('pole-zero plot')
H0=freqz(b,a,w);phi=unwrap(angle(H0));
Kn=real((10^(dB/20)-0.5*(H0+1))/(0.5*(1-H0)));
disp(['f0=' num2str(vf0(1,n)) ',k1=' num2str(k1n) ',k2=' num2str(k2n) ',K=' num2str(Kn)])
subplot(132),semilogx(faxis,phi),grid,axis([0 fs/2 -2*pi 0])
xlabel('f,Hz'),ylabel('phase,rad'),title('phase response')
Fn=0.5*(1+Kn+(1-Kn)*H0);
subplot(133),semilogx(faxis,20*log10(abs(Fn)+eps))
grid,axis([0 fs/2 -10 10])
xlabel('f,Hz'),ylabel('K(f), dB'),title('frequency response')
if n==1
F1=Fn;k11=k1n;k21=k2n;K1=Kn; str2=[s2 s4 s8];
elseif n==2
F2=Fn;k12=k1n;k22=k2n;K2=Kn;str2=[s2 s5 s8];
elseif n==3
F3=Fn;k13=k1n;k23=k2n;K3=Kn; str2=[s2 s6 s8];
elseif n==4
F4=Fn;k14=k1n;k24=k2n;K4=Kn; str2=[s2 s7 s8];
else
F5=Fn;k15=k1n;k25=k2n;K5=Kn;
end
end
F=F1.*F2.*F3.*F4.*F5;
figure(6)
semilogx(faxis,20*log10(abs(F)+eps))
grid,axis([0 fs/2 -10 10])
xlabel('f,Hz'),ylabel('K(f), dB'),title('total frequency response')
information5=str2mat('You can test the magnitide response equalizer',...
'with some real voice signal.',...
'Please,make your choice:',...
'1-Press 1 for testing the equalizer',...
'2-Press 2 for Exit.');
disp(information5)
choice3=input('Please,make your choice:');
switch choice3
case 1
load handel
sos=[0.5*((1+k21)+K1*(1-k21)) k11*(1+k21) 0.5*((1+k21)-K1*(1-k21)) 1 k11*(1+k21) k21
0.5*((1+k22)+K2*(1-k22)) k12*(1+k22) 0.5*((1+k22)-K2*(1-k22)) 1 k12*(1+k22) k22
0.5*((1+k23)+K3*(1-k23)) k13*(1+k23) 0.5*((1+k23)-K3*(1-k23)) 1 k13*(1+k23) k23
0.5*((1+k24)+K4*(1-k24)) k14*(1+k24) 0.5*((1+k24)-K4*(1-k24)) 1 k14*(1+k24) k24
0.5*((1+k25)+K5*(1-k25)) k15*(1+k25) 0.5*((1+k25)-K5*(1-k25)) 1 k15*(1+k25) k25];
figure(7)
subplot(121),psd(y);
title('Power Spectral Density Estimate via Wellch for the input signal')
y1=sosfilt(sos,y);
subplot(122),psd(y1);
title('Power Spectral Density Estimate via Wellch for the output signal')
information3=str2mat('You can hear the input and output signal.',...
'Please,make your choice:',...
'1-Press 1 for hearing the signals',...
'2-Press 2 for Exit.');
disp(information3)
choice2=input('Please,make your choice:');
switch choice2
case 1
wavplay(y)
wavplay(y1)
exitFlag=1;
case 2,exitFlag=1;
end
case 2,exitFlag=1;
end
case 3,exitFlag=1;
end
end
