%ECE 480 Final Project
%Steve Spiteri
%July 30, 2003

%%%%%%%%%%%%%%%%%%%%
%Anti-Aliasing LPF
%%%%%%%%%%%%%%%%%%%%

fsamp = 10e3;
fp = 4.98e3;
fs = 5.29e3;
ap = .5;
as = 40;
[n,wn] = ellipord(2*pi*fp, 2*pi*fs, ap, as, 's');
[num1, den1] = ellip(n, ap, as, wn, 's');
f = [0:100:10e3]';
h = freqs(num1, den1, 2*pi*f);
plot(f/(1e3), abs(h));
xlabel('Frequency in kHz'), ylabel('Amplitude'), title('Magnitude response in linear scale'), grid on;
figure;
att = -20*log10(abs(h));
plot(f/1e3, att);
xlabel('Frequency in kHz'), ylabel('Attenuation in dB'), title('Attenuation versus frequency');
axis([0, 10, 0, 50]), grid on;
figure;

%%%%%%%%%%%%%%%%%%%%
%Low Pass Filter #1
%%%%%%%%%%%%%%%%%%%%


%Filter Specifications
f1 = 262.5;      %Cutoff Frequency
fsamp = 1.2e4;  %Sampling Frequency

wc = f1 * 2 * pi/fsamp; %digital frequency
a1 = (1-sin(wc))/cos(wc);    %Find a
b = [1-a1] * [1,1]; %num vector of the trans
a = 2*[1,-a1]; %den vector of trans
f = [0:.1:5e3]';    %Frequency in Hz
omega = 2 * pi * f/fsamp;   %Digital Frequency
H = freqz(b,a,omega);   %Response
magdb = 20 * log10(abs(H));
plot(f,magdb), grid on
axis([0, 5000, -45, 1]), xlabel('Frequency in Hz'), ylabel('Magnitude Response in dB');
title('Magnitude Response of the Digital LPF #1');
figure;


%%%%%%%%%%%%%%%%%%%%
%Band Pass Filter #1
%%%%%%%%%%%%%%%%%%%%


%Filter Specifications
f1 = 262.5;   %Stopband
f2 = 675;   %Stopband
as = 40;    %Stopband attenuation in dB
flow = 362.5; %Lower passing frequency in Hz
fupper = 575; %Upper passing frequency in Hz
ap = 0.5;   %Passband attenuation in dB
fsamp = 1.2e4;    %Sampling frequency

%Specification in Digital Frequency
fomega = [f1, flow, fupper, f2] * 2 * pi/fsamp;  %Specification Vector

%Prewarping the digital filter
T = 2;  %Parameter for prewarp

%Analog filter specification
W_anal = 2/T * tan(fomega/2);

%Design of analog Elliptic filter
[N, wn] = ellipord(W_anal(2:3), W_anal(1:3:4), ap, as, 's');
%Note that 2N is the order (only for BPF and Notch filters)
[B,A] = ellip(N, ap, as, wn, 's');  %Analog filter coefficients

%Design of digital filter using the bilinear transformation
[b1,a1] = bilinear(B,A,1/T);
%b and a are the coefficients of the digital filter
f = [0:.1:5e3]';    %Frequency in Hz
omega = 2 * pi * f/fsamp;   %Digital Frequency
H = freqz(b1,a1,omega);   %Response
magdb = 20 * log10(abs(H));
plot(f,magdb), grid on
axis([200, 800, -45, 1]), xlabel('Frequency in Hz'), ylabel('Magnitude Response in dB');
title('Magnitude Response of the Digital Elliptic BPF #1');
figure;


%%%%%%%%%%%%%%%%%%%%
%Band Pass Filter #2
%%%%%%%%%%%%%%%%%%%%


%Filter Specifications
f1 = 575;   %Stopband
f2 = 1300;   %Stopband
as = 40;    %Stopband attenuation in dB
flow = 675; %Lower passing frequency in Hz
fupper = 1200; %Upper passing frequency in Hz
ap = 0.5;   %Passband attenuation in dB
fsamp = 1.2e4;    %Sampling frequency

%Specification in Digital Frequency
fomega = [f1, flow, fupper, f2] * 2 * pi/fsamp;  %Specification Vector

%Prewarping the digital filter
T = 2;  %Parameter for prewarp

%Analog filter specification
W_anal = 2/T * tan(fomega/2);

%Design of analog Elliptic filter
[N, wn] = ellipord(W_anal(2:3), W_anal(1:3:4), ap, as, 's');
%Note that 2N is the order (only for BPF and Notch filters)
[B,A] = ellip(N, ap, as, wn, 's');  %Analog filter coefficients

%Design of digital filter using the bilinear transformation
[b2,a2] = bilinear(B,A,1/T);
%b and a are the coefficients of the digital filter
f = [0:.1:5e3]';    %Frequency in Hz
omega = 2 * pi * f/fsamp;   %Digital Frequency
H = freqz(b2,a2,omega);   %Response
magdb = 20 * log10(abs(H));
plot(f,magdb), grid on
axis([500, 1400, -45, 1]), xlabel('Frequency in Hz'), ylabel('Magnitude Response in dB');
title('Magnitude Response of the Digital Elliptic BPF #2');
figure;


%%%%%%%%%%%%%%%%%%%%
%Band Pass Filter #3
%%%%%%%%%%%%%%%%%%%%


%Filter Specifications
f1 = 1200;   %Stopband
f2 = 2550;   %Stopband
as = 40;    %Stopband attenuation in dB
flow = 1300; %Lower passing frequency in Hz
fupper = 2450; %Upper passing frequency in Hz
ap = 0.5;   %Passband attenuation in dB
fsamp = 1.2e4;    %Sampling frequency

%Specification in Digital Frequency
fomega = [f1, flow, fupper, f2] * 2 * pi/fsamp;  %Specification Vector

%Prewarping the digital filter
T = 2;  %Parameter for prewarp

%Analog filter specification
W_anal = 2/T * tan(fomega/2);

%Design of analog Elliptic filter
[N, wn] = ellipord(W_anal(2:3), W_anal(1:3:4), ap, as, 's');
%Note that 2N is the order (only for BPF and Notch filters)
[B,A] = ellip(N, ap, as, wn, 's');  %Analog filter coefficients

%Design of digital filter using the bilinear transformation
[b3,a3] = bilinear(B,A,1/T);
%b and a are the coefficients of the digital filter
f = [0:.1:5e3]';    %Frequency in Hz
omega = 2 * pi * f/fsamp;   %Digital Frequency
H = freqz(b3,a3,omega);   %Response
magdb = 20 * log10(abs(H));
plot(f,magdb), grid on
axis([1100, 2750, -45, 1]), xlabel('Frequency in Hz'), ylabel('Magnitude Response in dB');
title('Magnitude Response of the Digital Elliptic BPF #3');
figure;


%%%%%%%%%%%%%%%%%%%%
%Band Pass Filter #4
%%%%%%%%%%%%%%%%%%%%


%Filter Specifications
f1 = 2450;   %Stopband
f2 = 5550;   %Stopband
as = 40;    %Stopband attenuation in dB
flow = 2550; %Lower passing frequency in Hz
fupper = 5450; %Upper passing frequency in Hz
ap = 0.5;   %Passband attenuation in dB
fsamp = 1.2e4;    %Sampling frequency

%Specification in Digital Frequency
fomega = [f1, flow, fupper, f2] * 2 * pi/fsamp;  %Specification Vector

%Prewarping the digital filter
T = 2;  %Parameter for prewarp

%Analog filter specification
W_anal = 2/T * tan(fomega/2);

%Design of analog Elliptic filter
[N, wn] = ellipord(W_anal(2:3), W_anal(1:3:4), ap, as, 's');
%Note that 2N is the order (only for BPF and Notch filters)
[B,A] = ellip(N, ap, as, wn, 's');  %Analog filter coefficients

%Design of digital filter using the bilinear transformation
[b4,a4] = bilinear(B,A,1/T);
%b and a are the coefficients of the digital filter
f = [0:.1:5e3]';    %Frequency in Hz
omega = 2 * pi * f/fsamp;   %Digital Frequency
H = freqz(b4,a4,omega);   %Response
magdb = 20 * log10(abs(H));
plot(f,magdb), grid on
axis([2300, 5700, -45, 1]), xlabel('Frequency in Hz'), ylabel('Magnitude Response in dB');
title('Magnitude Response of the Digital Elliptic BPF #4');

wavfile = 'hello';  %original wave file
[y,fs,wavbits] = wavread(wavfile);
sound(y,fs,wavbits);
plot(y);
title('Original Wave File');
xlabel('Samples, n'), ylabel('Amplitude');
Time = [0:length(y)-1]'*1/fs;
x = ([Time,y]);

%sound(output, fs);