%
%  EEE202Filter.m   Keith Holbert   November 2006
%
clear all
 
% create a time vector
power2 = 10;        % power of two
npts = 2^power2;    % number of data points per signal
period = 5;         % period of the signal (sec)
deltat = period / npts; % sampling interval (sec)
fsamp = 1/deltat;   % sampling frequency (Hz)
time = [0 : deltat : period-deltat];  % time vector (sec)

% Initialize a frequency vector (Hz) using linspace(a,b,n)
% which generates an n point vector from a to b inclusive.
nfreq = npts/2;
freq = linspace(0,1/(2*deltat),npts);

% Create a random input signal over the full time period
dc = 10;
for k = 1:npts
    sigin(k) = dc + rand(1);
end

%  Alternate signal of a sum of sines
%sigin = time * 0;
%for k = 1:nfreq
%    sigin = sigin + cos(2*pi*freq(k) * time);
%end
    
% plot the input signal for reference
subplot(2,1,1); 
plot(time,sigin);
title('Input Signal in Time Domain'); 
ylabel('Input Signal');
xlabel('Time (sec)');

% determine the input signal spectrum using the FFT
spectin = fft(sigin);

% plot the input signal spectrum for reference
subplot(2,1,2); 
loglog(freq(2:nfreq+1),abs(spectin(2:nfreq+1)));
title('Input Signal Spectrum'); 
ylabel('Input Magnitude');
xlabel('Frequency (Hz)');

figure;

% Now, let's use a simple low-pass filter to remove the
% high frequency components
% Set the cutoff frequency to one-eighth the sampling frequency
tau = 8/(2*pi*fsamp);
filter = (1/tau) ./ (j*2*pi.*freq + 1/tau);

% A Bode plot of the transfer function (log scale)
subplot(2,1,1); 
loglog(freq,abs(filter)); 
title('Filter Transfer Function Magnitude'); 
ylabel('Gain (units)');

subplot(2,1,2); 
semilogx(freq,180/pi*angle(filter));
title('Filter Transfer Function Phase'); 
ylabel('Angle (degrees)');

xlabel('Frequency (Hz)');
figure;

% The ouput signal is the product of the input signal spectrum
% and the filter transfer function in the frequency domain
spectout = spectin .* filter;

% Perform inverse FFT to obtain time domain signal
sigout = abs(ifft(spectout));

% Plots of the output signal and its spectrum
subplot(2,1,1); 
plot(time,sigout);
title('Output Signal in Time Domain'); 
ylabel('Output Signal');
xlabel('Time (sec)');

subplot(2,1,2); 
loglog(freq(2:nfreq+1),abs(spectout(2:nfreq+1)));
title('Output Signal Spectrum'); 
ylabel('Output Magnitude');
xlabel('Frequency (Hz)');