%INEL 4102 - BASIC ENGINEERING CIRCUITS
%CAUSAL FIRST ORDER RC-FILTER 
%Ideal Low-pass filter versus RC-filter
%Prof. Domingo Rodriguez
%
%PARAMETER SETTINGS**************************************
clear all                           %Clear all parameters
close all                           %Close all graphs and plots
Fs=5000;					     %Sampling Frequency
Ts=1/Fs;						  %Sampling Period or Sampling Time
V=0.040;						  %Time duration (in seconds) for each signal
R=1000;						     %Value of resistor
C=7.9577e-007;                %Value of the capacitor
fm=1/(2*pi*C*R);			 %Value of cut-off frequency of one-sided bandwidth
Wn=2*fm/Fs;					  %Normalized frequency
a=1/(R*C);					    %Time Constant Parameter
h0=a;							   %Initial Condition Parameter
t=0:Ts:V-Ts;				    %General time axis
N=length(t);                     %Number of samples
f=-Fs/2:Fs/N:Fs/2-Fs/N;	   %Frequency axis
%********************************************************
%RC-FILLTER 
hRC=h0*exp(-a*t);		%RC-Filter impulse response function
hmax=max(hRC);				%Maximum value of impulse response function
hRC=(hRC/hmax);			%Normalized impulse response function
fhRC=fft(hRC);				%Fourier Transform of impulse response function
sfhRC=fftshift(fhRC);	%Frequency shift for two sided spectrum plot
asfhRC=abs(sfhRC);		%Absolute value calculation
Hmax=max(asfhRC);			%Maximum value of frequency response function
asfhRC=(asfhRC/Hmax);	%Normilized frequency response function
%
%********************************************************
%MATLAB FIR1 FILTER
Wn=(2*fm)/Fs;				%Normilized cut-off frequency for ideal filter
h=fir1(N-1,Wn);			%Impulse response function of ideal filter
fh=fft(h);					%Fourier transform of impulse response function
sfh=fftshift(fh);			%Frequency shift for two sided spectrum plot
msfh=abs(sfh);				%Absolute value calculation
%
%********************************************************
%PLOTS
%********************************************************
plot(f,asfhRC,f,msfh,'r') 	%Plots of RC and Ideal filters
grid
xlabel('Frequency in Hz.')	%Horizontal axis
ylabel('Magnitude')			%Vertical axis
title('HL(f) and HRC(f)')	%Title of plot