%% Estimate the Power Spectrum Using dsp.SpectrumEstimator
% Alternately, you can compute the power spectrum of the signal using the |dsp.SpectrumEstimator| 
% System object. You can acquire the output of the spectrum estimator and store
% the data for further processing. To view other objects in the |Estimation| library, 
% type |help dsp| in the MATLAB&reg; command prompt, and click |Estimation|. 
%% Initialization
% Use the same source as in the previous section on using the  |dsp.SpectrumAnalyzer| 
% to estimate the power spectrum. The input sine wave has two frequencies: one 
% at 1000 Hz and the other at 5000 Hz. Initialize |dsp.SpectrumEstimator| to compute 
% the power spectrum of the signal using the filter bank approach. View the power 
% spectrum of the signal using the |dsp.ArrayPlot| object.

Fs = 44100;
Sineobject1 = dsp.SineWave('SamplesPerFrame',1024,'PhaseOffset',10,...
    'SampleRate',Fs,'Frequency',1000);
Sineobject2 = dsp.SineWave('SamplesPerFrame',1024,...
    'SampleRate',Fs,'Frequency',5000);

SpecEst = dsp.SpectrumEstimator('Method','Filter bank',...
    'PowerUnits','dBm','SampleRate',Fs,'FrequencyRange','onesided');
ArrPlot = dsp.ArrayPlot('PlotType','Line','ChannelNames',{'Power spectrum of the input'},...
    'YLimits',[-80 30],'XLabel','Number of samples per frame','YLabel',...
    'Power (dBm)','Title','One-sided power spectrum with respect to samples');
%% Estimation
% Stream in and estimate the power spectrum of the signal. Construct a |for|-loop 
% to run for 5000 iterations. In each iteration, stream in 1024 samples (one frame) 
% of each sine wave and compute the power spectrum of each frame. Add Gaussian 
% noise with mean at 0 and a standard deviation of 0.001 to the input signal.
%%
for Iter = 1:5000
    Sinewave1 = Sineobject1();
    Sinewave2 = Sineobject2();
    Input = Sinewave1 + Sinewave2;
    NoisyInput = Input + 0.001*randn(1024,1);
    PSoutput = SpecEst(NoisyInput);
    ArrPlot(PSoutput);
end

%% 
% Using the filter bank approach, the spectral estimate has a high resolution 
% and the peaks are precise with no spectral leakage.
%% Convert _x_-axis to Represent Frequency
% By default, the array plot shows the power spectral data with respect to the 
% number of samples per frame. The number of points on the x-axis equals the length 
% of the input frame. The spectrum analyzer plots the power spectral data with 
% respect to frequency. For a one-sided spectrum, the frequency varies in the 
% range [0 Fs/2]. For a two-sided spectrum, the frequency varies in the range 
% [-Fs/2 Fs/2]. To convert the x-axis of the array plot from sample-based to frequency-based, 
% do the following:
% 
% * Click on the *Configuration Properties* icon.
% * For a one-sided spectrum - On *Main* tab, set *Sample increment* to $Fs/FrameLength$ 
% and *X-offset* to 0.
% * For a two-sided spectrum - On *Main* tab, set *Sample increment* to $Fs/FrameLength$ 
% and *X-offset* to $-Fs/2$.
% 
% In this example, the spectrum is one-sided and hence, the *Sample increment* 
% and *X-offset* are set to 44100/1024 and 0, respectively. To specify the frequency 
% in kHz, set the *Sample increment* to 44.1/1024.
%%
ArrPlot.SampleIncrement = (Fs/1000)/1024;
ArrPlot.XLabel = 'Frequency (kHz)';
ArrPlot.Title = 'One-sided power spectrum with respect to frequency';

for Iter = 1:5000
    Sinewave1 = Sineobject1();
    Sinewave2 = Sineobject2();
    Input = Sinewave1 + Sinewave2;
    NoisyInput = Input + 0.001*randn(1024,1);
    PSoutput = SpecEst(NoisyInput);
    ArrPlot(PSoutput);
end


%% Live Processing
% The output of the |dsp.SpectrumEstimator| object contains the spectral data 
% and is available for further processing. The data can be processed in real-time 
% or it can be stored in the workspace.

%% 
% Copyright 2012 The MathWorks, Inc.