%480 Final Project
%Spring/Summer 2003
clear;
clc;
g1=1;   %gain 1
g2=1;   %gain 2
g3=1;   %gain 3
g4=1;   %gain 4
fsamp= 1e4;  %digital sampling
t=1/fsamp;  %digital period
T=2;         %peroid for warping and prewarping

%stop bands in format filter 1..4 frequency 1
f11=150;
f21=350;
f31=750;
f41=1550;

%stop bands in format filter 1..4 frequency 2
f12=450;
f22=850;
f32=1650;
f42=3250;

%-pass bands in format filter 1..4 frequency lo
f1lo=250;
f2lo=450;
f3lo=850;
f4lo=1650;

%-pass bands in format filter 1..4 frequency up
f1hi=350;
f2hi=750;
f3hi=1550;
f4hi=3150;

%att in dB
ap= 0.5;
as= 40;

%converting analog into dital
filt1=[f11, f1lo, f1hi, f12]*2*pi/fsamp;
filt2=[f21, f2lo, f2hi, f22]*2*pi/fsamp;
filt3=[f31, f3lo, f3hi, f32]*2*pi/fsamp;
filt4=[f41, f4lo, f4hi, f42]*2*pi/fsamp;

%prewarping
T=2;
filter_a1= 2/T*tan(filt1/2);
filter_a2= 2/T*tan(filt2/2);
filter_a3= 2/T*tan(filt3/2);
filter_a4= 2/T*tan(filt4/2);

%here come the filters
[N1, wn1] = ellipord (filter_a1(2:3), filter_a1(1:3:4),ap,as,'s');
[N2, wn2] = ellipord (filter_a2(2:3), filter_a2(1:3:4),ap,as,'s');
[N3, wn3] = ellipord (filter_a3(2:3), filter_a3(1:3:4),ap,as,'s');
[N4, wn4] = ellipord (filter_a4(2:3), filter_a4(1:3:4),ap,as,'s');
%note for bandpass fitlers

[B1,A1]=ellip(N1, ap, as, wn1, 's');
[B2,A2]=ellip(N2, ap, as, wn2, 's');
[B3,A3]=ellip(N3, ap, as, wn3, 's');
[B4,A4]=ellip(N4, ap, as, wn4, 's');

%frequency response of the four filters
freqs(B1,A1)
title('Frequency Response of 1st BPF');
figure(2)
freqs(B2,A2)
title('Frequency Response of 2nd BPF');
figure(3)
freqs(B3,A3)
title('Frequency Response of 3rd BPF');
figure(4)
freqs(B4,A4)
title('Frequency Response of 4th BPF');

figure(5)
%here is the filter conversion form analog to digital using the bilinear transform
[b1,a1]=bilinear(B1,A1,1/T); 
[b2,a2]=bilinear(B2,A2,1/T); 
[b3,a3]=bilinear(B3,A3,1/T); 
[b4,a4]=bilinear(B4,A4,1/T);
disp('Band pass filters complete')


%test to see if filters work correctly**************************************
disp('test start')

f=[0:.1:5e3]';
omega = 2*pi*f/fsamp;
H1=freqz(b1,a1,omega);
H2=freqz(b2,a2,omega);
H3=freqz(b3,a3,omega);
H4=freqz(b4,a4,omega);

dB_H1=20*log10(abs(H1));
dB_H2=20*log10(abs(H2));
dB_H3=20*log10(abs(H3));
dB_H4=20*log10(abs(H4));

subplot(2,2,1)
plot(f, dB_H1)
grid on;
axis([0,500,-45,1]);
xlabel('Frequency in Hertz');
ylabel('Magnitude in dB');
title('Mag response of digital BPF   1');

subplot(2,2,2)
plot(f, dB_H2)
grid on;
axis([0,900,-45,1]);
xlabel('Frequency in Hertz');
ylabel('Magnitude in dB');
title('Mag response of digital BPF   2');

subplot(2,2,3)
plot(f, dB_H3)
grid on;
axis([0,1700,-45,1]);
xlabel('Frequency in Hertz');
ylabel('Magnitude in dB');
title('Mag response of digital BPF   3');

subplot(2,2,4)
plot(f, dB_H4)
grid on;
axis([0,3300,-45,1]);
xlabel('Frequency in Hertz');
ylabel('Magnitude in dB');
title('Mag response of digital BPF   4');

figure(6)
subplot(2,2,1)
plot(f, -dB_H1)
grid on;
axis([0,500,-1,45]);
xlabel('Frequency in Hertz');
ylabel('Att in dB');
title('Att of digital BPF   1');

subplot(2,2,2)
plot(f, -dB_H2)
grid on;
axis([0,900,-1,45]);
xlabel('Frequency in Hertz');
ylabel('Att in dB');
title('Att of digital BPF   2');

subplot(2,2,3)
plot(f, -dB_H3)
grid on;
axis([0,1700,-1,45]);
xlabel('Frequency in Hertz');
ylabel('Att in dB');
title('Att of digital BPF   3');

subplot(2,2,4)
plot(f, -dB_H4)
grid on;
axis([0,3300,-1,45]);
xlabel('Frequency in Hertz');
ylabel('Att in dB');
title('Att of digital BPF    4');
figure(7)
disp('test complete')
%******************************************************************************
%continus to digital
fs = 4500;
fp = 3500;
[N,Wn] = ellipord(2*pi*fp, 2*pi*fs, ap, as, 's');
[NUM_aa, DEN_aa] = ellip(N, ap, as, Wn, 's');     %finding the transfor coeff
f=[0:1:8e3];
H= freqs(NUM_aa, DEN_aa, 2*pi*f);

subplot(3,1,1)
plot(f/1e3, abs(H));
xlabel ('Freq in kHz');
ylabel ('Amplitude');
title('Magnitude response in linear scale');
grid on;

att= -20*log10(abs(H)); %Att in Db
subplot(3,1,2)
plot(f/1e3 , att)
title ('Attenuation versus freq')
axis([0,10,0,50]);
grid on;

subplot(3,1,3)
plot(f/1e3,angle(H))
xlabel('Hz');
ylabel('Radians')
title('Phase response');
grid on;

figure
 
%***********************************************************************
%ad and da converter values
nb = 8;%**********************number of bits
vref= 5;
delta = (2^(1 -nb))*vref;
l1 = vref - delta;
l2 = -vref;
%wave**************************
wavfile = 'voice1';
[Y, FS, BITS]=wavread(wavfile);
plot(Y)
title ('Original wave data')
xlabel('nsamples')
ylabel('Amplitude')
figure
wavplay(Y, FS);

Time=[0:length(Y)-1]'*1/FS;
plot(Time,Y);
Data=[Time,Y];


figure
plot(OUTPUT(:,1),OUTPUT(:,2));
title ('Filtered wave data')
xlabel('nsamples')
ylabel('Amplitude')
wavplay(OUTPUT)


















