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

% The purpose ofthis file is to sample the sum of two previously filtered
% sinusoids at a rate of 44.1 kHz.  In order to do this, two sine waves will be
% created which have 44100 data points in them, with the number of cycles
% required for one second.  The sinusoids will then be composed of vectors that
% have been sampled at 44.1 kHz.  The next step is to quantize each magnitude
% of the sampled signal into a specific, discrete value.  Since each sample
% can be represented by a 16 bit binary number, only 65536 values are allowed.
% So, the magnitude will be normalized to the maximum signal value that can be
% represented. Then, the normalized value will be multiplied by 65535.  This
% number will be rounded to the nearest integer. This will be the value of the
% signal at that point.  This value will be stored in a vector describing the
% sampled sinusoid.  Each sinusoid will complete one full cycle every
% 44,100/sinusoid frequency samples 
% This is from the formula samples/cycle = samples/second * seconds/cycle





% first, we form a matrix for all the frequencies we want from 0 to 32000 Hz

w=0:20*2*pi:32000*2*pi;
numfreqs = length (w);

% Next, we build the inverse Chebyshev filter
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');

% finally, we create teh magnitude and phase shift vectors for each value of w
[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.

% sin3 is a 30kHz filtered sine wave.  It yields 30 thousand cylces every
% second. It repeats a cycle every 4.41/3 sample points, when sampling at 44.1
% kHz.  Notice that this will result in aliasing, since we are not obeying the
% Nyquist sampling rate requirement.  Therfore, the filter was designed to
% eliminate this signal as much as possible.

% shift3 is for the 30kHz signal
% shift2 is for the 10kHz signal
% shift1 is for the 1 kHz signal

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

shiftsin3 = .3*sin(2*pi*[0:44100]/(4.41/3) + shift3);
shiftsin2 = .1*sin(2*pi*[0:44100]/4.41 + shift2);
shiftsin1 = .1*sin(2*pi*[0:44100]/44.1 + shift1);

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


qfilt3 = round(filt3*65536)/65536;
qfilt2 = round(filt2*65536)/65536;
qfilt1 = round(filt1*65536)/65536;

x3 = qfilt1 + qfilt2 + qfilt3;




% remove after testing

plot (qfilt1(1:100));
ylabel ('amplitude of input 1kHz filtered signal in volts');
title(' filtered sine wave at 1kHz, sampled at 44.1kHz ');
xlabel ('sample points, time equal to (sample point#)/44,100');  
pause;

plot (qfilt2(1:100));
ylabel ('amplitude of input 10kHz filtered signal in volts');
title(' filtered sine wave at 10kHz sampled at 44.1kHz ');
xlabel ('sample points, time equal to (sample point#)/44,100');  
pause;
plot (qfilt3(1:100));
ylabel ('amplitude of input 30kHz filtered signal in volts');
title(' filtered sine wave at 30kHz sampled at 44.1kHz ');
xlabel ('sample points, time equal to (sample point#)/44,100');  
pause;

stem (x3(1:100));
ylabel ('amplitude of input filtered signal');
title('stem plot of three filtered sine waves at 1kHz, 10kHz, 30kHz  summed and sampled at 44.1kHz ');
xlabel ('sample points, time equal to (sample point#)/44,100');