% eda09_14
% example of Chebychev bandpass filter

% i used here as imaginary unit!
clear i;

% limits of bandpass filter
flow=5;
fhigh=10;

% standard time/frequency setup
N=1024;
Nf=N/2+1;
Dt=0.01;
tmax=Dt*(N-1); 
t=Dt*[0:N-1]'; 
fmax=1/(2.0*Dt);
Df=fmax/(N/2);
f=Df*[0:N/2,-N/2+1:-1]'; 
Dw=2*pi*Df;
w=Dw*[0:N/2,-N/2+1:-1]'; 

% input timeseries is a spike
din=zeros(N,1);
din(100)=1;

% filter
[dout, u, v] = chebyshevfilt( din, Dt, flow, fhigh);

% fft
dinspec = abs(fft(din));
doutspec = abs(fft(dout));

% plot input
figure(1);
clf;
subplot(2,2,1);
set(gca,'LineWidth',2);
hold on;
axis( [0.75, 1.25, -1, 1 ] );
plot( t, din, 'k-', 'LineWidth', 2);
xlabel('time, s');
ylabel('input');

% plot input spectrum
subplot(2,2,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, 1.1*max(dinspec) ] );
plot( f(1:Nf), dinspec(1:Nf), 'k-', 'LineWidth', 2);
xlabel('frequency, Hz');
ylabel('input spectrum');

% plot output
subplot(2,2,3);
set(gca,'LineWidth',2);
hold on;
axis( [0.75, 1.25, -1, 1 ] );
plot( t, dout, 'k-', 'LineWidth', 2);
xlabel('time, s');
ylabel('output');

% plot output spectrum
subplot(2,2,4);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, 1.1*max(dinspec) ] );
plot( f(1:Nf), doutspec(1:Nf), 'k-', 'LineWidth', 2);
xlabel('frequency, Hz');
ylabel('output spectrum');

% roots of filters
ru = roots(flipud(u));
Nu = length(u);
Nru = length(ru);
rv = roots(flipud(v));
Nv=length(v);
Nrv = length(rv);

% complex z-plane setup
Nz=100;
zmin=-2;
zmax=2;
dz=(zmax-zmin)/(Nz-1); 
zr=zmin+dz*[0:Nz-1]'; 
zi=zmin+dz*[0:Nz-1]';

% evaluate filter polynomials everywhere in complex z-plane
S = zeros(Nz,Nz);
for n = [1:Nz]
for m = [1:Nz]
    z = zr(n)+i*zi(m);
    uu=0.0;
    zpow=1;
    for k = [1:Nu]
        uu = uu + u(k)*zpow;
        zpow=zpow*z;
    end
    vv=0.0;
    zpow=1;
    for k = [1:Nv]
        vv = vv + v(k)*zpow;
        zpow=zpow*z;
    end
    SS = uu/vv;
    S(m,n)=abs(SS);
end
end   

% figure showing position of roots in the z-plane
% plot amplitude of polynomial
figure(2);
clf;
set(gca,'LineWidth',2);
hold on;
axis ij;
axis equal;
zmin=-2;
zmax=2;
axis([zmin, zmax, zmin, zmax]);
bw=0.9*(256-linspace(0,255,256)')/256;
colormap([bw,bw,bw]);
range=max(max(S))-min(min(S));
if( range==0 )
    range=1;
end
imagesc( [zmin, zmax], [zmin, zmax], (S-min(min(S)))/range);
xlabel('real z');
ylabel('imag z');


% plot origin
plot(0,0,'ok','LineWidth',2);
z=exp(-i*w*Dt);
plot(real(z),imag(z),'k-','LineWidth',2);

% plot zeros and poles
for k = [1:Nru]
    plot(real(ru(k)),imag(ru(k)),'k*','LineWidth',2);
end
for k = [1:Nrv]
    plot(real(rv(k)),imag(rv(k)),'k+','LineWidth',2);
end

