function PSM=PlotSegFreqSpect(y,NumSegments,MaxHarmonic,logmode)
% function PSM=PlotSegFreqSpect(y,NumSegments,MaxHarmonic,logmode)
% Plots the "Short-time Fourier transform (STFT)" of y. Breaks y into
% 'NumSegments' equal-length segments, multiplies each by a Hanning window,
% computes the power spectrum of each segment, and plots the result of the
% first 'MaxHarmonic' Fourier components as a contour plot. If the number
% of segments and of data points is such that the last segment is
% incomplete, it is discarded. Returns the power spectrum matrix
% (time-frequency-amplitude) as a matrix PSM of size NumSegments X
% MaxHarmonic. If logmode=1, it computes and plots the base10 log of the
% amplitudes. Requires the "hamming.m" function from the signal processing toolkit.
% Version 1.3 (c) Tom O'Haver January, 2020.
%
% EXAMPLE 1: Segmented power spectrum of passing automobile horn
% See https://terpconnect.umd.edu/~toh/spectrum/CaseStudies.html#Harmonic
% load horngoby
% PSM=PlotSegFreqSpect(doppler,40,350,0); % 40 segments, 240 harmonics
% figure(2);mesh(PSM); % Figure window 2 shows the 3-D mesh plot of the
% % power spectrum matrix, which can be rotated by dragging the pointer).
%
% EXAMPLE 2: Segmented power spectrum of spoken phrase "Testing 1, 2, 3"
% load testing123; % Download from https://terpconnect.umd.edu/~toh/spectrum/
% t=0:1/44001:1.5825; 
% y=testing123(:,1); % Select one of the two stereo channels.
% PSM=PlotSegFreqSpect(y,40,160,0); % 40 segments, 160 harmonics.
%
% EXAMPLE 3: Two Gaussians with constant sinewave noise, displayed in log
% mode.
% x=1:12000;
% y=90*gaussian(x,3000,500)+90*gaussian(x,7500,500);
% y=y+sin(x/2); % Add sinewave "noise"
% clf
% PSM=PlotSegFreqSpect(y,40,30,1);
%
% EXAMPLE 4: Sinusoidal "chirp", log mode
% 
% clf;
% x=1:.1:2000;
% f=.00002.*x;
% z=sin(2*pi.*f.*x);
% subplot(2,1,1);
% plot(x,z);
% subplot(2,1,2);
% NumSegments=7;
% MaxHarmonic=25;
% PSM=PlotSegFreqSpect(z,NumSegments,MaxHarmonic,1);
%
% EXAMPLE 5: Morse code S-O-S at fixed frequency
% load sos % Download sos.mat from https://terpconnect.umd.edu/~toh/spectrum/
% PlotSegFreqSpect(sos,100,200,0);
%
% EXAMPLE 6: Detecting an invisible signal
% Detection a weak signal totally buried in a combination of random and
% periodic noise. The simulated signal consists of three components: a
% pair of weak Gaussian peaks, a strong variable-frequency sine wave, and
% random white noise. The Gaussian peaks are invisible in the raw signal
% but show as yellow blobs at x= 2 and 4, y=1, in the segmented spectrum,
% clearly distinct from the much stronger diagonal stripe of yellow from
% the swept sine wave.
% User-adjustable parameters
% h=0.15; % Height of the Gaussian peak
% p1=500;  % x-axis position of the Gaussian peak 1
% p2=1000;  % x-axis position of the Gaussian peak 2
% w=250; % width of the Gaussian peak
% noise=0.3;
% smoothwidth=1100;
% NumSegments=7;
% MaxHarmonic=15;
% x=1:.1:2000;
% f=.00002.*x;
% z=sin(2*pi.*(max(f)-f./2).*x); % Compute swept sine wave
% g=h.*gaussian(x,p1,w)+h.*gaussian(x,p2,w);
% y=g+z+noise.*randn(size(g));
% figure(1);
% clf;plot(x,y);
% title('Raw signal')
% xlabel('time')
% figure(2);
% clf
% PSM=PlotSegFreqSpect(y,NumSegments,MaxHarmonic,1);
%
% Related functions: The built-in Matlab function S = stft(X) returns the
% short-time Fourier transform (stft) of the vector or matrix X.
% FouFilter.m, PlotSegFrequencySpectrum.m,SegmentedFouFilter.m, isignal.m. 

sizey=size(y);if sizey(1)>sizey(2);y=y';end
if NumSegments<2;disp('Number of segments must be >2');NumSegments=2;end
if MaxHarmonic<1;disp('Number of harmonics must be >1');end
PSM=[];
ly=length(y);
SegLength=round(ly./NumSegments);
fy=zeros(1,ly);
sy=zeros(1,ly);
for Segment=1:NumSegments
    startindex=(1+(Segment-1)*SegLength);
    endindix=startindex+SegLength-1;
    if endindix>ly,endindix=ly;end
    indexrange=startindex:endindix;
    fy(indexrange)=fft(y(indexrange).*hamming(length(indexrange))');
    sy(indexrange)=fy(indexrange) .* conj(fy(indexrange));
    try
      PSM=[PSM;sy(indexrange)];
    catch 
      NumSegments=NumSegments-1;    
    end
end
if logmode
   PSM=log10(PSM);
end

contour(PSM(1:NumSegments,1:MaxHarmonic)',200); 
% contourf(PSM(1:NumSegments,1:MaxHarmonic)',20); 
xlabel('Time Segments')
ylabel('Frequency (Harmonics)')
if logmode
   title('Segmented Power Spectrum, log z-axis amplitude scale')
else
   title('Segmented Power Spectrum, linear z-axis amplitude scale')
end
grid
