echo off;clear;clc;close all;warning off all
%%%%%% Example from Matlab
Fs = 1000;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector
gTt = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); %%% The sum of two sinosoidal signal
%%%% One positive rectangular pulse
gTt = [zeros(1,round(length(t)/2)) ones(1,round(length(t)/10)) zeros(1,round(8*length(t)/20))];
%%%% Two rectangular pulses with zero DC (zero mean)
gTt = [zeros(1,round(length(t)/2)) ones(1,round(length(t)/20)) -ones(1,round(length(t)/20)) zeros(1,round(8*length(t)/20))];
% gTt = [zeros(1,round(8*length(t)/20)) -ones(1,round(length(t)/20)) ones(1,round(length(t)/20)) -ones(1,round(length(t)/20)) ....
%       ones(1,round(length(t)/20)) zeros(1,round(8*length(t)/20))];
figure(1)
plot(1000*t(1:end),gTt(1:end),'linewidth',3)
title('Input signal');xlabel('\itt \rm(milliseconds)');ylabel('\itg_T \rm( \itt\rm )')
set(gcf,'Color',[1 1 1]);set(gca,'FontSize',14);
Gf = fft(gTt);
Gf2 = abs(Gf/L); %%% Note the L scaling
Gf1 = Gf2(1:L/2+1);
f = Fs*(0:(L/2))/L;ff = Fs*(-L/2:L/2 - 1)/L; %%%% f is one sided frequency range, ff is two sided
figure(2)
plot(f,Gf1,'linewidth',2);set(gcf,'Color',[1 1 1]);set(gca,'FontSize',14); 
title('Single-sided amplitude spectrum of \itG_T \rm( \itf\rm )')
xlabel('\itf \rm(Hz)');ylabel('|\itG_T \rm( \itf\rm )|')
figure(3)
plot(ff,ifftshift(Gf2),'linewidth',2)  %%% ifftshift is used to correctly position the spectrum 
set(gcf,'Color',[1 1 1]);set(gca,'FontSize',14);
title('Double-sided amplitude spectrum of \itG_T \rm( \itf\rm )')
xlabel('\itf \rm(Hz)');ylabel('|\itG_T \rm( \itf\rm )|')
%%%% The effect of filtering
fc = 100; Cf = [ones(1,fc) zeros(1,length(f) - fc)];%Cf = 1;%%% Produces abnormal results
sigma = 40;%Cf = exp(-f.^5/sigma^2); %%% Exponential low pass filter
fc = 1000;Cf = 1./(1 + j*f.^2/fc);Cf = ones(1,length(f));
figure(4)
semilogx(f,Cf,'linewidth',2);set(gcf,'Color',[1 1 1]);set(gca,'FontSize',14); 
title('Bandlimited channel response');xlabel('\itf \rm(Hz)');ylabel('|\itC\rm ( \itf \rm)|')

Cf = Cf(1:end-1); % set correct length L/2
Gff = Gf.*[Cf fliplr(Cf)]; % double sided filter

figure(5)
%plot(ff,Sff) % abs + ifftshift required
plot(ff,abs(ifftshift(Gff))/L,'linewidth',2);set(gcf,'Color',[1 1 1]);set(gca,'FontSize',14);
title('Signal spectrum after passing through the channel');xlabel('\itf \rm(Hz)');ylabel('|\itG_T \rm( \itf\rm )|')
gtf = ifft(Gff);

figure(6)
%plot(abs(stf(1:end))) % abs not required
plot(1000*t,gtf,'-b','linewidth',3);hold on;plot(1000*t(1:end),gTt(1:end),'--r','linewidth',3);hold off
set(gcf,'Color',[1 1 1]);set(gca,'FontSize',14);
title('Input signal and signal after passing through bandlimited channel');xlabel('\itt \rm(milliseconds)')
ylabel('\itg_f \rm( \itt \it)')
%%% Energy of transmitted signal 
Et = 1e-3*gTt*gTt'
Er = 1e-3*gtf*gtf'
