%% Estimate the Power Spectrum Using dsp.SpectrumAnalyzer
% To view the power spectrum of a signal, you can use the |dsp.SpectrumAnalyzer|
% System object(TM).
% You can change the dynamics of the input signal and see the effect
% those changes have on the power spectrum of the signal in real time.

%% Initialization
% Initialize the sine wave source to generate the sine wave and the
% spectrum analyzer to show the power spectrum of the signal.
% The input sine wave has two frequencies: one at 1000 Hz and the other at
% 5000 Hz. Create two |dsp.SineWave| objects, one to generate the 1000 Hz
% sine wave and the other to generate the 5000 Hz sine wave.
Fs = 44100;
Sineobject1 = dsp.SineWave('SamplesPerFrame',1024,'PhaseOffset',10,...
'SampleRate',Fs,'Frequency',1000);
Sineobject2 = dsp.SineWave('SamplesPerFrame',1024,...
'SampleRate',Fs,'Frequency',5000);
SA = dsp.SpectrumAnalyzer('SampleRate',Fs,'Method','Filter bank',...
'SpectrumType','Power','PlotAsTwoSidedSpectrum',false,...
'ChannelNames',{'Power spectrum of the input'},'YLimits',[-120 40],'ShowLegend',true);
%%
% The spectrum analyzer uses the filter bank approach to compute the power
% spectrum of the signal.
%% 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.
% To generate the input signal, add the two sine waves. The resultant
% signal is a sine wave with two frequencies: one at 1000 Hz and the
% other at 5000 Hz. Add Gaussian noise with zero mean and a standard
% deviation of 0.001. To acquire the spectral data for further processing,
% use the |isNewDataReady| and the |getSpectrumData| object functions.
% The variable |data| contains the spectral data that is displayed on the
% spectrum analyzer along with additional statistics about the spectrum.
data = [];
for Iter = 1:7000
Sinewave1 = Sineobject1();
Sinewave2 = Sineobject2();
Input = Sinewave1 + Sinewave2;
NoisyInput = Input + 0.001*randn(1024,1);
SA(NoisyInput);
if SA.isNewDataReady
data = [data;getSpectrumData(SA)];
end
end
release(SA);
%%
% In the spectrum analyzer output, you can see two distinct peaks: one at
% 1000 Hz and the other at 5000 Hz.
%%
% Resolution Bandwidth (RBW) is the minimum frequency bandwidth that
% can be resolved by the spectrum analyzer. By default, the
% |RBWSource| property of the |dsp.SpectrumAnalyzer| object is set to |Auto|.
% In this mode, RBW is the ratio of the frequency span to 1024. In a
% two-sided spectrum, this value is $\frac{F_{s}}{1024}$, while in a one-sided
% spectrum, it is $\frac{\frac{F_{s}}{2}}{1024}$. The spectrum analyzer in
% this example shows a one-sided spectrum. Hence, RBW is (44100/2)/1024 or
% 21.53Hz
%
% Using this value of $RBW$, the number of input samples required to compute
% one spectral update, $N_{samples}$ is given by the following equation:
% $N_{samples} = \frac{F_{s}}{RBW}$.
%
% In this example, $N_{samples}$ is 44100/21.53 or 2048 samples.
%
% $RBW$ calculated in the 'Auto' mode gives a good frequency resolution.
%
% To distinguish between two frequencies in the display, the distance
% between the two frequencies must be at least RBW.
% In this example, the distance between the two peaks is 4000 Hz, which is
% greater than $RBW$. Hence, you can see the peaks distinctly.
% Change the frequency of the second sine wave to 1015 Hz. The
% difference between the two frequencies is less than $RBW$.
release(Sineobject2);
Sineobject2.Frequency = 1015;
for Iter = 1:5000
Sinewave1 = Sineobject1();
Sinewave2 = Sineobject2();
Input = Sinewave1 + Sinewave2;
NoisyInput = Input + 0.001*randn(1024,1);
SA(NoisyInput);
end
release(SA);
%%
% The peaks are not distinguishable.
%
% <<../zoomedimage_15.png>>
%
%%
% To increase the frequency resolution, decrease $RBW$ to 1 Hz.
SA.RBWSource = 'property';
SA.RBW = 1;
for Iter = 1:5000
Sinewave1 = Sineobject1();
Sinewave2 = Sineobject2();
Input = Sinewave1 + Sinewave2;
NoisyInput = Input + 0.001*randn(1024,1);
SA(NoisyInput);
end
release(SA);
%%
% On zooming, the two peaks, which are 15 Hz apart, are now distinguishable.
%
% <<../zoomedimage_1.png>>
%%
% When you increase the frequency resolution, the time resolution decreases.
% To maintain a good balance between the frequency resolution and time
% resolution, change the |RBWSource| property to |Auto|.

%%
% During streaming, you can change the input properties or the spectrum
% analyzer properties and see the effect on the spectrum
% analyzer output immediately.
% For example, change the frequency of the second sine wave when the index
% of the loop is a multiple of 1000.
%
release(Sineobject2);
SA.RBWSource = 'Auto';
for Iter = 1:5000
Sinewave1 = Sineobject1();
if (mod(Iter,1000) == 0)
release(Sineobject2);
Sineobject2.Frequency = Iter;
Sinewave2 = Sineobject2();
else
Sinewave2 = Sineobject2();
end
Input = Sinewave1 + Sinewave2;
NoisyInput = Input + 0.001*randn(1024,1);
SA(NoisyInput);
end
release(SA);
%%
% While running the streaming loop, you can see that the peak of the second
% sine wave changes according to the iteration value. Similarly, you can
% change any of the spectrum analyzer properties while the simulation is
% running and see a corresponding change in the output.

%%
% Copyright 2012 The MathWorks, Inc.