%Set cutoff frequency and delay of the impulse response peak
Fc = 1;
tau = 5;

%Two-sided bandwidth of the lowpass filter
B = 2*Fc;

%Set up time array and sampling variables for the "ideal" lowpass filter
%Impulse response length >> 1/B to approximate the ideal filter
time1 = linspace(-50/B,50/B,4000)';
dt1 = time1(2) - time1(1);      %Sampling interval for ideal filter time array
Fs1 = 1/dt1;                    %Sampling frequency
dt2 = 1/(50*B);                 %Sampling interval for practical filter time array
time2 = (0:dt2:(2*tau))';       %Time array for the practical filter
Fs2 = 1/dt2;                    %Sampling frequency
N = length(time2);              %Number of samples in practical filter

%Compute impulse responses for ideal and practical filters
ht_ideal = B*sinc(B*time1);
ht_windowed = B*sinc(B*(time2 - tau)).*ones(N,1);       %with rectangular window
% ht_windowed = B*sinc(B*(time2 - tau)).*hanning(N);    %with hanning window

%I have evaluated the Fourier Transform using the FFT command. You will be
%making your own Fourier transform script, so this section does not apply
%to you
Nfft = 16384;
F_axis1 = -0.5*Fs1 + (0:(Nfft-1))*Fs1/Nfft;
F_axis2 = -0.5*Fs2 + (0:(Nfft-1))*Fs2/Nfft;
% Hf_ideal = fftshift(abs(fft(ht_ideal,Nfft)));
% Hf_ideal = Hf_ideal/max(Hf_ideal);
Hf_ideal = 1e-6*ones(Nfft,1);
Hf_ideal(abs(F_axis1) < Fc) = 1;
Hf_windowed = fftshift(abs(fft(ht_windowed,Nfft)));
Hf_windowed = Hf_windowed/max(Hf_windowed);

%Now plot the results
figure
subplot(2,1,1)
plot(time1,ht_ideal,'linewidth',2)
hold
plot([-tau -tau tau tau],[0 B B 0],'r:','linewidth',2)  %show the equivalent width of the practical filter
set(gca,'fontsize',14)
xlabel('Time (s)')
ylabel('h(t)')
subplot(2,1,2)
plot(F_axis1,20*log10(Hf_ideal),'linewidth',2)
axis([-4*Fc 4*Fc -50 5])
set(gca,'fontsize',14)
xlabel('Frequency (Hz)')
ylabel('|H(f)|')

figure
subplot(2,1,1)
plot(time2,ht_windowed,'linewidth',2)
set(gca,'fontsize',14)
xlabel('Time (s)')
ylabel('h(t)')
subplot(2,1,2)
plot(F_axis2,20*log10(Hf_windowed),'linewidth',2)
axis([-4*Fc 4*Fc -50 5])
set(gca,'fontsize',14)
xlabel('Frequency (Hz)')
ylabel('|H(f)|')










