7.5.m

把各阶导函数的ESD作为a的函数

%
% FUNCTION 7.5:"cp0702_Gaussian_derivatives_peak
% frequency"
%
% Analysis of peak frequency of the first 15 derivatives of
% the Gaussian pulse as a funtion of the shape factor
%
% 'smp' samples of the Gaussian pulse are considered in
% the time interval 'Tmax-Tmin'
%
% The function receives as input;
% 1) the minimum value of the shape factor 'alphamin'
% 2) the increase step 'alphastep'
% 3) the number of values to be investigated
%    'N_alphavalues'
%
% The function computes the ESDs of
% the first 15 derivatives of the Gaussian pulse for the
% 'alpha' value received as input,and then evaluates and
% plots the peak frequency for each derivative
%
% Programmed by Luca De Nardis

funvtion cp0702_Gaussian_derivatives_peak_frequency...
    (alphamin,alphastep,N_alphavalues)

% --------------------------
% Step Zero-Input parameters
% --------------------------

smp=1024;             % number of samples
alpha=alphamin;       % Gaussian pulse form factor
Tmin=-4e-9;           % lower time limit
Tmax=4e-9;            % upper time limit

dt=(Tmax-Tmin)/smp;   % sampling period
fs=1/dt;              % sampling frequency
N=smp;                % number of samples(i.e.,size of
                      % the FFT)
df=1/(N*dt);          % fundamental frequency

x=linspace(Tmin,Tmax,smp);
F=figure(1);

for j=1:N_alphavalues
    
    factor(j)=alpha;
    for i=1:15;
       
% ----------------------------------------------
% Step One - Amplitude-normalized pulse waveform
%            in the time domain
% ----------------------------------------------

    derivative(i,:)=...
        cp0702_analytical_waveforms(x,i,alpha);
    detivative(i,:)=derivative(i,:)/...
        max(abs(detivative(i,:)));
    
% ----------------------------------------------------
% Step Two - Analysis in the frequency domain and peak
%            frequency evaluation
% ----------------------------------------------------

        % double-sided MATLAB amplitude spectrum
        X=fft(derivative(i,:),N);
        % conversion from MATLAB spectrum to Fourier
        % spectrum
        x=X/N;
        % double-sided ESD
        E = fftshift(abs(X).^2/(df^2));
        % single-sided ESD
        Ess = 2* E((N/2+1):N);
        
        positivefrequency=linspace(0,(fs/2),N/2);
        
        % evaluation of the peak frequency (frequency at
        % which the ESD assumes the maximum value)
        [peak,peakelementindex]=max(Ess);
        
   peakfrequency(i,j) = ...
      positivefrequency(peakelementindex);
   end
    % increase of alpha value for the next step
    alpha = alpha + alphastep;
end

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

PT=plot(factor,peakfrequency');
set(PT,'LineWidth',[2]);
AX=gca;
set(AX,'FontSize',12);
X=xlabel('\alpha [s]');
set(X,'FontSize',14);
Y=ylavel('Peak frequency [Hz]');
set(Y,'FontSize',14);
axis([2e-10 14e-10 0 12e9]);
derivebehaviour = {'Increasing differentiation order'};
text(5e-10,6e9,derivebehaviour,'BackgroundColor',...
    [1 1 1]);
text(1.15e-9,0.5e9,'1^{st} derivative',...
    'BackgroundColor',[1 1 1]);
text(1.15e-9,2.4e9,'15^{th} derivative','BackgroundColor',
[1 1 1]);