%% This code generates a square wave with a given frequency and duty cycle
% and then filters the signal with a first order RC low pass filter
%%
clear;
clf;
    
fo=10;        % Square wave fundamental frequency in Hz
To=1/fo;
T_disp=5*To;  % duration of the waveform to be displayed
N=2^16;       % length of the time array
Ts=T_disp/(N-1);  % sampling interval  
t=0:Ts:T_disp;    % time array for 5 cycles of the square wave

%% Generate and Plot a square wave.
duty=50;      % duty cycle is the percent of the period where the signal
              % is positive
x=square(2*pi*fo*t,duty); % generate a square wave with frequency fo and 
                          % a given duty cycle  
figure(1)
subplot(3,1,1);
plot(t,x);
ylim([ -1.2 1.2]);

title('Square Wave, Time Domain')
xlabel('Time, seconds')

%% Compute the fft of the square wave.

X_fft=fft(x, N);           % compute fft 
% The +ve frequency components of the signal are represented by the first
% half of the array(n=1 to N/2) , 
% The -ve frequency component(N/2+1 to N) of the signal are represented 
% by the second half
% fftshift(X) rearranges the outputs of fft,moving the -ve frequency 
% components to the beginning of the array
% It is useful for visualizing a Fourier transform with the zero-frequency 
% component in the middle of the spectrum.

X1_fft=fftshift(X_fft);

% generate the frequency axis
fs=1/Ts;            % sampling frequency
df=fs/(N);          % frequency spacing

Freq=(-fs/2):df:fs/2-df; % frequency array
X1_max=max(abs(X1_fft));
X1_norm=X1_fft/X1_max;   % normalized spectrum
subplot(3,1,2);
plot(Freq,abs(X1_norm));
title('Square wave, frequency domain')
xlabel('Frequency(Hz)')
xlim([-10*fo 10*fo]);

%% Transform to time domain
%Now transform back to the time domain and plot to demonstrate that the
%fft function followed by the ifft function gives the original time
%domain waveform .
x2=ifft(X_fft);
subplot(3,1,3);
plot(t,x2)
xlabel('Time, seconds')
title('Square wave, time domain, using ifft')
ylim([ -1.2 1.2]);

%% Filter transfer Function 
%Now filter the frequency domain version with a low-pass filter.
%H(w)=-Rf/Rs*/(1+jw*Rf*C)
%Assume that R=100k and C=0.1uF, 
% Note that the filtering is performed on the frequency domain version of 
% the signal and considering the positive and the negative frequency 
% portions of the signal. The spectrum X1_fft is used which is centered
% around zero frequency 

figure(2)
Rf=1e5   ;
Rs=1e5;
C=1e-7;
% filter transfer function H(w)
Hw=-Rf/Rs*1./(1+1j*2*pi*Freq*Rf*C); % Freq is the frequency axis with positive 
                                    % and negative frequencies
subplot(3,1,1)
plot(Freq, abs(Hw))
title('Magnitude of Filter Transfer Function H(w)')
xlabel('Frequency (Hz)')
xlim([-10*fo 10*fo])

%% Filtering the input signal
% calculating the fft of the filtered square wave in the frequency domain
X1_fft_filtered=X1_fft.*Hw; % fft of the filtered signal
X1_filtered_norm=X1_fft_filtered/X1_max; % normalizing filtered spectrum
subplot(3,1,2)
plot(Freq,abs(X1_filtered_norm))
title(' Spectrum of Filtered Signal')
xlabel('Frequency (Hz)')
xlim([-10*fo 10*fo ])

%% 
%% Converting the filtered signal back to the time domain 
% in order to find the inverse fft, the spectrum of the signal has to be
% modified so that the zero frequency is at the beginning of the array, the
% command used is ifftshift which moves the first N/2 elements of the array
% to the end of the array. It swaps the right and left portion of the array
x_filt=ifft(ifftshift(X1_fft_filtered)); % filtered signal, time domain
subplot(3,1,3);
x_filt=real(x_filt);
plot(t,x_filt)
title('Output Filtered Waveform in the Time Domain')
xlabel('Time(s)')


