fprintf('\n\n =============================================================== \n');
fprintf(' semestralni prace CFS, Petr Sladek CTU FEE, sladep1: \n');
fprintf(' ------------------------------------------------------\n');
fprintf(' BAND-PASS DIGITAL FILTER:\n');
fprintf(' f0 = 770 Hz \n BWf0 = +/-8 Hz, passband \n BWf0 = +/-60 Hz, stopband \n');
fprintf(' ripple attenuation ap = 0.5dB, passband\n attenuation = 30dB, stopband \n');
fprintf(' approximation: Cauer/elliptic\n fsamp = 3kHz\n ');
fsamp = 3000;
f0 = 770;
fl = (f0-8);
fr = (f0+8);
BW = 8+8;
stopatt = 30;
passripple = 0.5;
flstop = f0-60;
frstop = f0+60;
flDstop = f0-100;
frDstop = f0+100;
typfs = 'elliptic';
wl = 2*pi*fl;
wr = 2*pi*fr;
wlstop = 2*pi*flstop;
wrstop = 2*pi*frstop;
f = flDstop:0.5:frDstop;
w = 2*pi*f;
if ((strcmp(typfs,'elliptic'))|(strcmp(typfs,'cheby1')))
else
fprintf('Wrong filter approximation specified!');
quit cancel;
end
if strcmp(typfs,'elliptic')
[Nfilt Wn] = ellipord([wl wr],[wlstop wrstop], passripple-0.02,stopatt, 's');
[Zs Ps Ks] = ellipap(Nfilt, passripple, stopatt);
end
if strcmp(typfs,'cheby1')
[Nfilt Wn] = cheb1ord([wl wr],[wlstop wrstop], passripple-0.02,stopatt, 's');
[Zs Ps Ks] = cheb1ap(Nfilt, passripple);
end
[Bs As Bz Az] = proto2dbp(Zs,Ps,Ks,f0,BW,fsamp);
[Bd Ad] = ellip(Nfilt,passripple,stopatt,[wl/(pi*fsamp) wr/(pi*fsamp)]);
fprintf('\n Filter designed, order %d:, approximation: %s \n',Nfilt,typfs)
Bz
Az
[SOS sosG] = tf2sos(Bz,Az,'up','two');
sosB(:,1:3) = SOS(:,1:3);
sosA(:,1:3) = SOS(:,4:6);
fprintf('and it''s second-order sections (SOS) form sosB (numerator), sosA(denom.), and overall gain sosG: \n')
sosB
sosA
sosG
[sosN foo] = size(SOS);
sosH=ones(sosN,length(w));
figure(2); clf;
for i=1:sosN,
subplot(sosN,2,2*(i-1)+1);
sosH(i,:)=freqz(sosB(i,:),sosA(i,:),w/fsamp);
plot(f,20*log10(abs(sosH(i,:))),'r'); grid on;
title(sprintf('second-order sections (SOS) magnitude of block %d, overall gain: %0.5g dB ',i,20*log10(sosG)));
xlabel('frequency [Hz]');
ylabel('magnitude [dB]');
subplot(sosN,2,2*i);
plot(f,angle(sosH(i,:)),'b'); grid on;
title(sprintf('second-order sections (SOS) phase of block %d',i));
xlabel('frequency [Hz]');
ylabel('phase [rad]');
end
Hsosall = prod(sosH)*sosG;
fprintf(' Evaluating impulse response, takes some time. Be patient please. \n');
[hs ts]=impulse(tf(Bs,As));
[hz tz]=impulse(tf(Bz,Az,1/fsamp));
fprintf(' DONE.\n');
figure(1); clf;
Hs=freqs(Bs,As,w);
subplot(2,1,1);
plot(f,20*log10(abs(Hs)),'r:'); grid on; hold on; ylim([-50 0]);
subplot(2,1,2);
plot(f,angle(Hs),'r:'); grid on; hold on;
Hz=freqz(Bz,Az,w/(fsamp));
subplot(2,1,1);
plot(f,20*log10(abs(Hz)),'b');
subplot(2,1,2);
plot(f,angle(Hz),'b');
Hd = freqz(Bd,Ad,w/(fsamp));
subplot(2,1,1);
plot(f,20*log10(abs(Hd)),'g--');
subplot(2,1,2);
plot(f,angle(Hd),'g--');
subplot(2,1,1);
plot(f,20*log10(abs(Hsosall)),'m-.');
subplot(2,1,2);
plot(f,angle(Hsosall),'m-.');
subplot(2,1,1);
legend('analog filter (r)','digital -manual impulse inv. (b)','digital -automated dsn (g) ellip only','SOS form (m)');
title(sprintf('magnitude of analog+digital filter: %s, order: %d, sampling f: %d Hz',typfs,Nfilt,fsamp));
xlabel('frequency [Hz]');
ylabel('magnitude [dB]');
plot([flDstop flstop flstop],[-stopatt -stopatt 0],'k','LineWidth',4);
plot([frDstop frstop frstop],[-stopatt -stopatt 0],'k','LineWidth',4);
plot([fl fl fr fr],[-stopatt -passripple -passripple -stopatt],'k','LineWidth',2);
plot([flstop frstop],[0 0],'k','LineWidth',2);
subplot(2,1,2);
legend('analog filter (r)','digital -manual impulse inv. (b)','digital -automated dsn (g) ellip only','SOS form (m)');
title(sprintf('phase of analog+digital filter: %s, order: %d, sampling f: %d Hz',typfs,Nfilt,fsamp));
xlabel('frequency [Hz]');
ylabel('phase [rad]');
figure(4); clf;
Hz=freqz(Bz,Az,w/(fsamp));
plot(f,20*log10(abs(Hz)),'b'); grid on; hold on; ylim([-(10*passripple) 0]); xlim([flstop frstop]);
plot([flDstop flstop flstop],[-stopatt -stopatt 0],'k','LineWidth',4);
plot([frDstop frstop frstop],[-stopatt -stopatt 0],'k','LineWidth',4);
plot([fl fl fr fr],[-stopatt -passripple -passripple -stopatt],'k','LineWidth',2);
title('detail view of digital filter - magnitude');
xlabel('frequency [Hz]');
ylabel('magnitude [dB]');
figure(3); clf;
subplot(2,1,1);
plot(ts,hs,'r'); grid on;
title('impulse response of analog denormalized prototype');
ylabel('magnitude');
xlabel('time [s]');
subplot(2,1,2);
stem(tz,hz.*fsamp,'b'); grid on;
title(sprintf('impulse response of digital designed filter, order: %d, fsampling: %d',Nfilt,fsamp));
xlabel('time [s]');
ylabel('magnitude');
figure(2);
fprintf('\n ====\n design&preview completed. Answers: sladek@asix.cz\n ==== \n');
===============================================================
semestralni prace CFS, Petr Sladek CTU FEE, sladep1:
------------------------------------------------------
BAND-PASS DIGITAL FILTER:
f0 = 770 Hz
BWf0 = +/-8 Hz, passband
BWf0 = +/-60 Hz, stopband
ripple attenuation ap = 0.5dB, passband
attenuation = 30dB, stopband
approximation: Cauer/elliptic
fsamp = 3kHz
Filter designed, order 2:, approximation: elliptic
Bz =
0.0301 0.0055 0.0590 0.0055 0.0302
Az =
1.0000 0.1657 1.9593 0.1619 0.9541
and it's second-order sections (SOS) form sosB (numerator), sosA(denom.), and overall gain sosG:
sosB =
0.2396 0.0768 0.2395
0.5857 -0.0807 0.5865
sosA =
1.0000 0.1166 0.9766
1.0000 0.0491 0.9770
sosG =
0.2149
Evaluating impulse response, takes some time. Be patient please.
DONE.
====
design&preview completed. Answers: sladek@asix.cz
====