%% Streaming Power Spectrum Estimation Using Welch's Method
% Compute the power spectrum estimate of a time-domain input signal using 
% the <docid:dsp_ref#buzhyoj-1 Spectrum Estimator> block. 
% The block uses one of the following methods to compute the power spectrum
% estimate:
%
% * Welch's method of averaged modified periodograms
% * Filter bank method
%
% This example uses the Welch's method of averaged modified periodograms.
% For an example that uses the filter bank-based spectrum estimation method, see 
% <docid:dsp_ug#example-HighResolutionSpectralAnalysisExample High
% Resolution Spectral Analysis>. The same example also shows the comparison
% between the filter bank estimator and the Welch-based spectral estimator.
% Generally, filter bank-based spectrum estimation yields better resolution 
% with less spectral leakage, more accurate peaks and a more accurate noise floor.
%
% For more details on algorithm of these two methods, see the 'Algorithms' 
% section in the <docid:dsp_ref#buzhyoj-1 Spectrum Estimator> block.
% 
% The Spectrum Estimator block is useful if you need direct access to the 
% estimated spectrum (rather than just visualize it). The output power 
% spectrum may be used as an input to other blocks in your model, or may be 
% logged to the workspace for post-processing. To visualize the spectra,
% use the <docid:dsp_ref#btqjkji Spectrum Analyzer> scope block.
%
% Copyright 2020 The MathWorks, Inc.
%% Welch's Method of Averaged Modified Periodograms
% In the Welch method, the input time-domain data is partitioned into data 
% segments based on the selected window length and overlap percentage. 
% This stage is implemented using a <docid:dsp_ref#f3-907916 Buffer> block.
% A window is applied to each segment, and then an averaged periodogram is 
% computed based on the windowed sequences. This stage is implemented using
% a <docid:dsp_ref#bt0cv78-1 dsp.SpectrumEstimator> System object(TM).
% The length of the data segments and the choice of the window determine the 
% estimate's resolution bandwidth (RBW), which is the smallest positive frequency 
% that can be resolved in the power spectral estimate.
%% Specifying Window Length
% The |dspstreamingwelch| model shown below uses a Welch Spectrum Estimator 
% block to estimate the spectrum of a noisy chirp signal sampled at 44100 Hz. 
% The power spectrum estimate is displayed using an <docid:dsp_ref#buzorp5-1 Array Plot> scope. 
% The peak value of the spectrum, as well as the frequency at which the peak occurs, 
% are detected and displayed on the scope. The estimate's RBW is also displayed. 
% Moreover, a Spectrum Analyzer scope block is included for comparison and 
% validation purposes.
open_system('dspstreamingwelch')

%%
% The block's frequency resolution method is set to |Window length|. The 
% window length is set to 1024. The FFT length, _NFFT_ , is equal to the 
% window length. The data is windowed using a Chebyshev window with a 
% sidelobe attenuation of 60 dB. The frequency range is one-sided. In this 
% case, the length of the spectrum estimate is $NFFT/2+1 = 513$ and is computed 
% over the interval [0 Hz,22050 Hz]. The *Sample increment* property of the 
% Array Plot scope is accordingly set to $Fs/NFFT = 44100/(1024 * 1000)$ , 
% where the increment is divided by 1000 to scale the frequency units to kHz. 
% You can access the scope's *Sample increment* property by opening its 
% Configuration properties window.
%
% The resolution bandwidth is given by:
% 
% $$RBW = enbw(chebwin(N,SL)) * Fs / N$$ 
%
% where _N_ is the window length, _enbw_ is a function that computes the 
% window's equivalent noise bandwidth, _SL_ is the sidelobe attenuation of 
% the selected Chebyshev window and _Fs_ is the sample rate. In this case, 
% RBW is equal to 65.38 Hz.
%%
% When you simulate the model, you can verify that the displayed RBW value 
% is equal to the one shown on the lower bar of the Spectrum Analyzer scope. 
% Moreover, the two blocks give the same peak measurements.
sim('dspstreamingwelch')
%% Specifying Non-Zero Overlap
% The model in the previous section had zero-overlap. In the |dspstreamingwelch_overlap| model, 
% we use a Welch Estimation block with an overlap of 50%. Since other model 
% parameters are identical to the previous section, the RBW is unchanged and 
% is equal to 65.38 Hz.
% With a window length of 1024 and an overlap percentage of 50%, 512 
% input samples are required to form a new data segment. Since the input 
% data is of length 1024, each new data frame yields two new periodograms, 
% and the block's output port runs at a rate twice as fast as the input port.
% 
% Note that the Welch estimate block does not have zero latency in this case. 
% The first spectrum estimate output is based on the buffer's initial condition, 
% which is equal to eps. In order to match the spectrum and measurements of 
% the Spectrum Analyzer scope, we therefore insert a delay block at the 
% input of the Spectrum Analyzer.
%
% The results of the Spectrum Analyzer and Welch estimate block may be validated 
% by simulating the model.
sim('dspstreamingwelch_overlap')
%% Specifying RBW
% In the |dspstreamingwelch_rbw| model, the *Frequency Resolution Method*
% parameter is set to |RBW|. RBW source is Auto. In this mode, similar to the 
% Spectrum Analyzer scope block, the resolution bandwidth is chosen such 
% that there are 1024 RBW intervals over the specified Frequency Span. 
% Since the span in this case is 22050 HZ, the RBW is 21.53 Hz.
%
% The window length used to buffer the data is iteratively computed to yield 
% the desired RBW. The window length in this case is equal to 3073. To verify this value, 
% we can compute the RBW that results from this window length:
%
% $$ RBW = enbw(hann(3073)) * 44100 / 3073 = 21.53 Hz $$
%
% Note that a Hann window is used in this model. In this case, the FFT length, NFFT, 
% is odd and equal to 3073 (the window length). Since the frequency range is 
% one-sided, the spectrum estimate is of length (NFFT + 1)/2 and is computed 
% over the interval [0,44100/2). The *Sample increment* property of the Array Plot 
% scope is set to $Fs/NFFT = 44100/(3073 * 1000)$ KHz.
%
% Again, the results of the Spectrum Analyzer and Welch estimate block can 
% be validated by simulating the model.
sim('dspstreamingwelch_rbw')
