% Digital Elliptic Band-pass design:
% Petr Sladek, @ CTU FEE

% semestralni prace CFS
% obsah:
% navrhnout jak matlab fcemi tak "rucne" z H(s) filtr s paramatery:
% impulsni invariance
%
% 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

% note: unfortunately there is big problem with using cells with publish() -
%       multiple plot() calls cause stupid multiple output... :-(.


% print input techdata:
%
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;

% analog. prototype
% analog prototype defs:
f0 = 770;             %center freq /f0\
fl = (f0-8);          %left /f0 pass
fr = (f0+8);          %right f0\ pass
BW = 8+8;             %bandwidth
stopatt = 30;         %stopband attenuation dB
passripple = 0.5;     %dB passband ripple
flstop = f0-60;       %stopband "left"  freq.
frstop = f0+60;       %stopband "right"  freq.
flDstop = f0-100;     %draw limits
frDstop = f0+100;
typfs = 'elliptic';    %specify approximation, {elliptic, cheby1}


% =====================================================================

% omegas:
wl = 2*pi*fl;
wr = 2*pi*fr;
wlstop = 2*pi*flstop;
wrstop = 2*pi*frstop;
f = flDstop:0.5:frDstop; %plot purposes freq_min->freq_max
w = 2*pi*f;


% evaluate min order of filter:
% passripple-0.02 - correction to match requirements
if ((strcmp(typfs,'elliptic'))|(strcmp(typfs,'cheby1')))
else
  fprintf('Wrong filter approximation specified!');
  quit cancel;
end

% analog prototype:
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

% manual impulse invariance (see fcn proto2dbp):
% returns analog denormalized filter also.
[Bs As Bz Az] = proto2dbp(Zs,Ps,Ks,f0,BW,fsamp);

% Matlab's automated digital bandpass design:
[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


% convert B, A into second-order sections (SOS) form:
% Bz/Az is input:
[SOS sosG] = tf2sos(Bz,Az,'up','two');

% extract B,A,sosG from SOS:
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

% SOS freq-response
[sosN foo] = size(SOS);
sosH=ones(sosN,length(w)); %malloc
figure(2); clf;

% view SOS freq response:
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;
    % sosH can be modified:
    % why  "*power(sosG,1/sosN)" ??? - because sosG is overall gain
    % H = sosG*H1*H2*...Hnm, power(sosG,1/sosN) splits sosG equally.
    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


% overall SOS form freq response product(sosH_i):
Hsosall = prod(sosH)*sosG;


% evaluate impulse response:
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');

% draw analog, digital (manual impulse invariance), digittal (automated)


% draw analog filter
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;

% draw digital filter (manual impulse inv)
%figure(1); problem with publish
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');




% draw direct discrete filter (automated design):

%figure(1); problem with publish
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--');

%figure(1); problem with publish
subplot(2,1,1);
plot(f,20*log10(abs(Hsosall)),'m-.');
subplot(2,1,2);
plot(f,angle(Hsosall),'m-.');

% show legend, title and so on:
%figure(1);
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]');
% draw filter specification:
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]');

% draw digital filter - detail view
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]);
% some mist to detail view..
% figure(4); problem with publish
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]');


% draw impulse response

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); %recall to move at the end of html> publish



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
 ==== 

Published with MATLAB® 7.0