% Reconstruction of a Gaussian peak from its lowest frequency components
% 
clear
format compact
maxx=1; % maximum value of x
dx=.0005; % sampling interval (.02 to .0001 or less)
x=0:dx:maxx; % Create independent variable x
N=length(x);
NumberOfPoints=N;
SamplingInterval=dx;
SamplingRate=1/dx;
NyquistFrequency=1/(2*dx);
width=max(x)/20; % Chanage the denominator to change the width inversely.
Noise=0; % Standard deviation of random noise (try 0.1 to 10)
g=gaussian(x,max(x)/4,width)+2*gaussian(x,max(x)/2,width)+3*gaussian(x,max(x)/1.4,width);
y=g+Noise.*randn(size(x)); % Add them up and add optional noise

% Upper panel - Gaussian signal peak
subplot(3,1,1);
plot(x,y,'linewidth',2);
grid
title('Multiple Gaussians signal')
xlabel(['Time, sec.     Center peakwidth (FWHM)= ' num2str(width)   '     Noise:'  num2str(Noise)  '     Total area: '  num2str(trapz(x,y))])
ylabel('Amplitude')
grid

% Compute Fourier spectrum "fs"
fy=fft(y);
sy=sqrt(fy .* conj(fy));
plotrange=1:length(fy)/2;
realsy=real(sy(plotrange));
f=((plotrange-1)./range(x)); % Frequency axis
grid

% lower panel - Reconstructed signal
flimit=1.25/width;  % Highest frequency to include in reconstruction
if flimit>floor(N/2),flimit=floor(N/2);end
for fmax=1:flimit 
fy=fft(y);
lfft=length(fy); % Length of the FFT
fy(fmax:lfft-fmax)=0; % All frequencies between fmax and lfft-fmax in the fft are set to zero.
ify=real(ifft(fy)); % Real part of the inverse fft
subplot(3,1,3);
plot(x,ify,'b','linewidth',2)
title(['Signal sequentially reconstructed from its '  num2str(fmax)  ' lowest frequencies.'])
xlabel(['time, sec.  Center peakwidth = ' num2str(halfwidth(x,ify))  '    Total area: '  num2str(trapz(x,ify)) ])
ylabel('amplitude')
grid

% Middle panel -  Frequency spectrum 
subplot(3,1,2);
plot(f,abs(realsy),[fmax fmax],[0 max(abs(realsy))],'linewidth',2)
axis([0 50 0 max(realsy) ])
grid
title('Frequency spectrum of the signal')
xlabel('Frequency, Hz. Signal dominates the lower frequencies; noise domimates the higher frequencies')
ylabel('amplitude')
drawnow
pause(.5) % Number of seconds pause between steps
end

function g = gaussian(x,pos,wid)
%  gaussian(x,pos,wid) = gaussian peak centered on pos, half-width=wid
%  x may be scalar, vector, or matrix, pos and wid both scalar
%  T. C. O'Haver, 1988
% Examples: gaussian([0 1 2],1,2) gives result [0.5000    1.0000    0.5000]
% plot(gaussian([1:100],50,20)) displays gaussian band centered at 50 with width 20.
g = exp(-((x-pos)./(0.60056120439323.*wid)) .^2);
end

function [FWHM,slope1,slope2,hwhm1,hwhm2]=halfwidth(x,y,xo)
% function [FWHM,slope1,slope2,x1,x2]=halfwidth(x,y,xo) computes the full
% width at half maximum of the maximum of any shape peak that peaks near xo
% and has a zero baseline. If xo is omitted, it computes the halfwidth from
% the maximum y value. Not highly accurate if the function is too sparsely
% sampled. Also optionally returns the leading and trailing edge slopes
% (slope1 and slope 2, respectively) and the leading and trailing edge
% half-width-at-half-maximum, if those output arguments are incuded. 
% Tom O'Haver (toh@umd.edu) Version 2, July 2016
%
% Example 1:
% x=-5:.01:5;
% y=sinc(x).^2;
% FWHM=halfwidth(x,y,0) % Area of central peak
% FWHM=halfwidth(x,y,1.5) % Area of positive side peak
%
% Example 2:
% x=[0:.1:10];
% W=3; % W is the true half-width
% y=gaussian(x,5,W);
% FWHM=halfwidth(x,y,5)
% FWHM=halfwidth(x,y) % gives the same result since there is only one peak
%
% Example 3: % Returns NaN because signal does not contain the halfwidth
% x=[0:.1:1];
% y=exp(-(x).^2)
% halfwidth(x,y) 
% 
% Example 4: Also returns the leading and trailing edge slopes (+1 and -1)
% and the leading and trailing edge half-width-at-half-maximum
% x=[0:.1:10];
% y=triangle(x,5,1);
% [FWHM,slope1,slope2,x1,x2]=halfwidth(x,y)
%
% Example 5. Comparison of Gaussian (g) and BiGaussian peaks (bg)
% x=[1:100];
% g=gaussian(x,50,20);
% bg=bigaussian(x,50,20,3);
% plot(x,g,x,bg)
% [GaussianFWHM,slope1,slope2,hwhm1,hwhm2]=halfwidth(x,g)
% [BiGaussianFWHM,slope1,slope2,hwhm1,hwhm2]=halfwidth(x,bg)

if nargin==2,xo=x(val2ind(y,max(y)));end

try   
    indmax=val2ind(x,xo);
    maxy=y(indmax);
    oy=y-maxy/2;
  
    n=indmax;
    while oy(n)>0
        n=n-1;
    end
    x1=interp1([y(n) y(n+1)],[x(n) x(n+1)],maxy/2);
    slope1=(y(n+1)-y(n))./(x(n+1)-x(n));
    
    n=indmax;
    while oy(n)>0
        n=n+1;
    end
    x2= interp1([y(n-1) y(n)],[x(n-1) x(n)],maxy/2);
    slope2=(y(n-1)-y(n))./(x(n-1)-x(n));
    
    FWHM=x2-x1;
    hwhm1=xo-x1;
    hwhm2=x2-xo;
catch me
    FWHM=NaN;
    slope1=NaN;
    slope2=NaN;
    hwhm1=NaN;
    hwhm2=NaN;
    me=me
end
end

function [index,closestval]=val2ind(x,val)
% Returns the index and the value of the element of vector x that is closest to val
% If more than one element is equally close, returns vectors of indicies and values
% Tom O'Haver (toh@umd.edu) October 2006
% Examples: If x=[1 2 4 3 5 9 6 4 5 3 1], then val2ind(x,6)=7 and val2ind(x,5.1)=[5 9]
% [indices values]=val2ind(x,3.3) returns indices = [4 10] and values = [3 3]
dif=abs(x-val);
index=find((dif-min(dif))==0);
closestval=x(index);
end