%
% FUNCTION 1.2 : "cp0101_bandwidth"
%
% Evaluates the bandwidth of the input 'signal' with sampling period
% 'dt'
% Bandwidth is evaluated according to the given 'threshold' (in dB)
% 'BW' is the bandwidth
% 'f_high' is the higher limit
% 'f_low' is the lower limit
%
% Programmed by Guerino Giancola
%
width = 1e-1;      % width of the rectangle [s]
points = 1000;     % number of samples for representing the rectangle  
f0 = 1e3;            % carrier frequency [Hz] 
dt = width / points;       % sampling period
signal = zeros(1,5*points);
signal(2*points:3*points-1)=ones(1,points);
mod=cos(2.*pi.*f0.*linspace(1,5*width,5*points));
signal=signal.*mod;
time=linspace(-2*width,3*width,5*points);

function[ss_E,f_high,f_low,BW]=...
    cp0101_bandwidth(signal,dt,threshold);
% ---------------------------------------------
% Step One - Evaluation of the single-sided ESD
% ---------------------------------------------

        % sampling frequency
N = length(signal);  % number of samples (i.e., size of the FFT)
T = N * dt;          % time window
df = 1 / T;          % fundamental frequency

X = fft(signal); % double-sided MATLAB amplitude spectrum
X = X/N;         % conversion from MATLAB spectrum to fourier spectrum
ds_E = abs(X).^2/(df^2);      % double-sided ESD
ss_E = 2.*ds_E(1:floor(N/2)); % single-sided ESD

% ------------------------------------------------
% Step Two - Evaluation of the frequency bandwidth
% ------------------------------------------------

[Epeak,index] = max(ss_E);      % Epeak is the peak value of the ESD
f_peak = index * df;            % peak frequency   

Eth = Epeak*10^(threshold/10);  % Eth is the value of the ESD
                                % corresponding to the given
                                % threshold

% iterative algorithm for evaluating high and low frequencies

imax = index;
E0h = ss_E(index);

while (E0h>Eth)&(imax<=(N/2))

    imax = imax + 1;
    E0h = ss_E(imax);

end % while E0h > Eth

f_high = (imax-1) * df;             % high frequency

imin = index;
E0l = ss_E(index);

while (E0l>Eth)&(imin>1)&(index>1)

    imin = imin - 1;
    E0l = ss_E(imin);

end % while E0l > Eth

f_low = (min(index,imin)-1) * df;   % low frequency

% end of iterative algorithm

BW = f_high - f_low;            % signal frequency bandwidth

fprintf('\nFrequency Bandwidth = %f [Hz]\nHigh Frequency = %f [Hz]\nLow Frequency =  %f [Hz]\n',BW,f_high,f_low);

% -----------------------------
% Step Three - Graphical output
% -----------------------------

figure(2)

frequency=linspace(0,fs/2,length(ss_E));
PF=plot(frequency,ss_E);
set(PF,'LineWidth',[2]);
L1=line([f_high f_high],[min(ss_E) max(ss_E)]);
set(L1,'Color',[0 0 0],'LineStyle',':')
L1=line([f_low f_low],[min(ss_E) max(ss_E)]);
set(L1,'Color',[0 0 0],'LineStyle',':')
L1=line([f_low f_high],[Eth Eth]);
set(L1,'LineWidth',[2],'Color','red','LineStyle',':')
axis([0.8*f_low 1.2*f_high -0.1*Epeak 1.2*Epeak]);
AX = gca;
set(AX,'FontSize',12);
T=title('Frequency domain');
set(T,'FontSize',14);
X=xlabel('Frequency [Hz]');
set(X,'FontSize',14);
Y=ylabel('Single-Sided ESD  [V^2s/Hz]');
set(Y,'FontSize',14);