%% Amplitude and Power Spectra of Two Separate Signals
% The amplitude and power spectra of nonstationary signals often provide no indication of how the frequency content of the signal changes with time. This frequency variatiom often contains critical information about the signal and process studied in applications. The magnitude spectrum, therefor, cannot suitably represent all the information contained in a nonstationary signal.
% Similarly, the Fourier transform is also inadequate for nonstationary signals, because their spectral representation requires that the variables of time and frequency be mutually exclusive. The evaluation of one value of the spectrum S(f) at frequency f requires the complete knowledge of the time domain signal s(t)  over -inf to inf.
% This limitation of the classical representations are well illustrated by the fact that we can find totally different signals related to different physical phenomena, such as a linear FM signal and a limited impulse, which have similar magnitude spectra.

N = 1200;
Fs = 700;
t = (0:N-1)/Fs;
rng default
x_a = 30*chirp(t,0,1,Fs/10);
SNR = 20;
noise = randn(size(x_a))*std(x_a)/db2mag(SNR);
x_a = x_a + noise;
plot(t,x_a);
title('Chirp Signal');
xlabel('Time (s)');
ylabel('Amplitude');

% aa

periodogram(x_a,[],[],Fs);
title('Periodogram Power Spectrum Estimate of a Noisy Linear Chirp')

t = 0:1/Fs:2;
x_b = 500*sinc(250*(t-1));
SNR = 18;
noise = randn(size(x_b))*std(x_b)/db2mag(SNR);
x_b = x_b + 2.5*noise;
plot(t,x_b)
grid on

%

figure;
periodogram(x_b,[],[],Fs);
title('Periodogram Power Spectrum Estimate of a Noisy Sinc')

%
%
%% Mesh plots of noisy chirp and noisy sinc in the time frequency domain
% Both the noisy linear frequency modulated chirp and the noisy sinc have very similar magnitude spectra.  To be able to differentiate between them, you would need to know the phases of these two signals. The sinc function is zero phase, whereas the frequency of the chirp is a function of time. You can reveal the time and frequency dependence of a signal by plotting a joint time frequency representation using a time frequency distribution. Time frequency representations also reveal whether the signal is monocomponent or multicomponent, a fact that cannot be revealed by conventional spectral analysis, especially when individual components are also time-varying.
% It is evident that nonstationary signals are more adequately represented through time frequency distributions. Time frequency representations are a natural extension of both time and frequency domain processing, using representations of signals in a complete space that can display all the signal information received from one sensor in a more accessible and interpretable way. Such representations are intended to provide a distribution of signal energy simultaneously over both time and frequency domains.

[p_a,f_a,t_a] = pspectrum(x_a,Fs,'spectrogram');
[p_b,f_b,t_b] = pspectrum(x_b,Fs,'spectrogram');
mesh(f_a',t_a',20*log10(abs(p_a))')
ylabel('Time (s)')
xlabel('Frequency (Hz)')
zlabel('Power (dB)')
zlim([-30 60])
colormap hsv

% aa

figure;
mesh(f_b',t_b',20*log10(abs(p_b))')
ylabel('Time (s)')
xlabel('Frequency (Hz)')
zlabel('Power (dB)')
zlim([-30 70])
colormap hsv

%% Time Frequency view of transient signals
% Visualize an interference narrowband signal embedded within a broadband signal.
% Generate a chirp sampled at 1 kHz for 500 seconds. The frequency of the chirp increases from 180 Hz to 220 Hz during the measurement.

Fs = 1000;
t = (0:1/Fs:500)';
x = chirp(t,180,t(end),220) + 0.15*randn(size(t));

% The signal also contains a 210 Hz sinusoid. The sinusoid has an amplitude of 0.05 and is present only for 1/6 of the total signal duration.

idx = floor(length(x)/6);
x(1:idx) = x(1:idx) + 0.05*cos(2*pi*t(1:idx)*210);

% Compute the spectrogram of the signal. Restrict the frequency range from 100 Hz to 290 Hz. Specify a time resolution of 1 second. Both signal components are visible.

pspectrum(x,Fs,'spectrogram', ...
    'FrequencyLimits',[100 290],'TimeResolution',1)
colormap default

% Compute the power spectrum of the signal. The weak sinusoid is obscured by the chirp.

pspectrum(x,Fs,'FrequencyLimits',[100 290])

% Compute the persistence spectrum of the signal. Now both signal components are clearly visible.

pspectrum(x,Fs,'persistence', ...
    'FrequencyLimits',[100 290],'TimeResolution',1)

%
% You can also visualize the two signal components using a constant-q stationary Gabor transform as shown below.

cqt(x,'SamplingFrequency',Fs,'FrequencyLimits',[100 290],'Window','hamming','BinsPerOctave',96)

% Copyright 2019 The MathWorks, Inc.