% Kurt Sanner
% Elec 301 Group Project
% December 15, 1999
% infilter.m

%timestart is the first time value for the x-axis
% maxtime is the largest time value for the x-axis
timestart = 0;
maxtime = 31999;

% sin1 repeats every 1000 units of time.  It yields 32 full cycles on the
% graph.  sin1 corresponds to a 30kHz sine wave.

sin1 = .3*sin(2*pi*[timestart:maxtime]/1000);

% sin2 repeats every 3000 units of time.  It yields 10 2/3 full cylces on the
% graph.  sin2 corresponds to a 10kHz sine wave.
sin2 = .1*sin(2*pi*[timestart:maxtime]/3000);

% sin3 repeats every 30000 units of time. It yields 1 1/15 full cycle on the
% graph. sin3 corresponds to a 1 kHz sine wave.
sin3 = .1*sin(2*pi*[timestart:maxtime]/30000);

% x1 is the unfiltered input sum of the three sinusoids.
x1 = sin1 + sin2 + sin3;

plot(x1);
ylabel ('amplitude of signal before filtering');
title(' sum of three sine waves at 1kHz, 10kHz, 30kHz, simulating an analog signal');
xlabel ('time, measured in microseconds');  


pause;

% x2 is the ideal result of the filtered input, which eliminates the high
% frequency component.
x2 = sin2 + sin3;

plot(x2)
ylabel ('amplitude of combined signals');
title(' plot of two lower frequency sine waves, simulating a perfectly filtered signal');
xlabel ('time, measured in microseconds') ; 

hold on;
pause;


% This part  creates an inverse Chebyshev filter
% The filter is a lowpass filter with the following properties:
%    Maximum decibel decrease at the passband is 2 decibels
%    The passband passband 2 decibel decrese is set at 17kHz
%    Minimum decibel decrease in the stopband is 50dB
%    The stopband 50 decibel decrease is set at 22kHz
%
% The inverse Chebyshev filter that results will be a eigth order filter.
% Although this filter will not be cheap, it is assumed to be used
% in a recording studio for optimal performance requirements.
format long g;
w=0:20*2*pi:32000*2*pi;
numfreqs = length (w);
Gp = 2;
Gs = 30;
Wp = 2*pi*18000;
Ws = 2*pi*22000;
[n, Wc] = cheb2ord (Wp, Ws, Gp, Gs, 's');
[num, den] = cheby2 (n, Gs, Wc, 's');
[MAG, PHASE] = bode(num, den, w);

%x3 is the filtered input

% there should be 1601 individual frequencies w for which the magnitude and
% phase shift have been calculated.  We need the particular magnitude and
% phase shift for the frequency of the sinusoid that is part of the input.

shift1 = 2*pi/360 * PHASE (round(numfreqs * 30000/32000));
shift2 = 2*pi/360 * PHASE (round(numfreqs * 10000/32000));
shift3 = 2*pi/360 * PHASE (round(numfreqs * 1000/32000));

shiftsin1 = .3*sin(2*pi*[0:31999]/1000 + shift1);
shiftsin2 = .1*sin(2*pi*[0:31999]/3000 + shift2);
shiftsin3 = .1*sin(2*pi*[0:31999]/30000 + shift3);

filt1 = MAG (round(numfreqs * 30000/32000)) * shiftsin1;
filt2 = MAG (round(numfreqs * 10000/32000)) * shiftsin2;
filt3 = MAG (round(numfreqs * 1000/32000)) * shiftsin3;
x3 = filt1 + filt2 + filt3;

plot (x3, 'r');

ylabel ('amplitude of signal after filtering');
title(' comparison of signal filtered with an inverse Chebyshev Filter(red) vs. ideal (blue)');
xlabel ('time, measured in microseconds');                    
hold off;