%% Measure Frequency Response by Using Spectrum Analyzer
% 
% This example shows how to emulate a spectrum analyzer that measures the
% frequency response of a continuous-time system driven by a complex
% sinusoidal signal. In this Simulink(R) model, the Plant block describes a
% second-order resonant system of the form 
%
% $$\frac{\omega^2}{s^2 + 2\zeta\omega s + \omega^2}$$
%
% with a natural frequency of $\omega = 2 \pi \cdot 150$ (or 150 Hz) and a
% damping ratio of $\zeta = 0.3$. Because the system is underdamped ($\zeta
% < 1$), the system has two complex conjugate poles. In a real-world
% application, replace this block with a subsystem composed of xPC D/A and
% A/D blocks that measure the response of a device under test (DUT).
%
% The spectrum analyzer consists of these components:
% 
% * The Sinusoidal Generator block, which produces a complex sinusoidal
% signal of increasing frequency. Inside the block, a Stateflow(R) chart
% uses temporal logic to iterate over a range of frequencies.
% * The Stateflow chart Analyzer, which calculates the frequency response
% (magnitude and phase angle) of the system at a specified frequency. The
% chart registers changes in the frequency by using change detection logic.
% * The Stateflow chart Unwrap, which processes the measured phase angle
% and unwraps the result so that there are no sharp jumps between $\pi$ and
% $-\pi$.
%
% The spectrum analyzer displays the measured frequency response as a pair
% of discrete Bode plots.

% Copyright 2008-2020 The MathWorks, Inc.

%% 

% Set up charts with white background
set_param(0,'ExportBackgroundColorMode','white');

mdl = 'sf_spectrum_analyzer';
open_system(mdl);
close_system([mdl '/Measured Frequency Response']);

%% Generate Sinusoid Signal
% The Sinusoid Generator block has two outputs:
% 
% * A scalar |f| that represents the current frequency.
% * A complex signal |y| that has a frequency of |f|.
% 
% You can control the behavior of the block by modifying its mask dialog
% box parameters. For example, the default values specify a sinusoidal
% signal with frequencies between 10Hz and 1000Hz. The block holds each
% frequency value for 0.2 seconds, and then increments the frequency by a
% step of 10Hz.
%
% <<../sf_spectrum_analyzer_signal_gen_mask.png>>
% 
% To control the timing of the signal generator, the block contains a
% Stateflow chart that applies absolute-time temporal logic.

open_system([mdl '/Sinusoid Generator/signal generator']);

%%
% During simulation, the chart goes through these stages:
%
% * *Initialize Frequency*: The default transition sets the signal
% frequency |f| to the value of the parameter |fstart|. Specify the value
% of |fstart| in the Sinusoid Generator mask dialog box.
% * *Generate Signal*: While the |Generate| state is active, the chart
% produces the complex signal |y = exp(2*pi*f*t*1i)| based on the frequency
% |f| and the time |t| since the last change in frequency. To determine the
% elapsed time (in seconds) since the state became active, the chart calls
% the temporal logic operator |<docid:stateflow_ref#temporalCount
% temporalCount>|.
% * *Update Frequency*: After the |Generate| state is active for |delay|
% seconds, the chart transitions out of the state, increases the frequency
% |f| by |fstep|, and returns to the |Generate| state. To determine the
% timing of the frequency updates, the chart calls the temporal logic
% operator |<docid:stateflow_ref#after after>| and checks the sign of the
% imaginary part of the signal before and after the current time step. The
% sign checks ensure that the output signal completes a full cycle before
% the frequency is updated, preventing large changes in the output signal.
% Specify the values of |fstep| and |delay| in the Sinusoid Generator mask
% dialog box.
% * *Stop Simulation*: When the frequency |f| reaches the value of the
% parameter |fstop|, the chart transitions to the |Stopped| state. The
% active state output |Stopped| triggers a <docid:simulink_ref#f8-894249>
% block and the simulation ends. Specify the value of |fstop| in the
% Sinusoid Generator mask dialog box.

%% Calculate Frequency Response
% The Analyzer chart takes the complex signal |u| and the output |yp| from
% the Plant block and computes the magnitude and phase of the plant output.

open_system([mdl '/Analyzer']);

%%
% For each frequency value, the chart maintains two running sums:
%
% * The sum |y| approximates the integral of the product of the plant
% output |yp| (a real number) and the complex signal |u|.
% * The sum |yn| approximates the integral of the product of the plant
% input |imag(u)| and the complex signal |u|. This integral represents the
% accumulation of a hypothetical plant with a unit transfer function.
%
% To detect a change in frequency, the chart uses the operator
% |<docid:stateflow_ref#hasChanged hasChanged>| to guard the
% self-transition on the state |Analyze|. When the frequency changes, the
% actions on this transition compute the magnitude and phase of the plant
% output by normalizing |y| with respect to |yn|.

%% Unwrap Measured Phase Angle
% The Unwrap chart prevents the measured phase angle from changing by more
% than $\pi$ radians in a single time step.

open_system([mdl '/Unwrap']);

%%
% In this chart, the transitions test the change in the input |u|
% before computing the new output value |y|.
%
% * If the input increases by more than $\pi$ radians, the chart offsets the
% output by $-2\pi$ radians.
% * If the input decreases by more than $-\pi$ radians, the chart
% offsets the output by $+2\pi$ radians.

%% Examine Simulation Results
% When you simulate the model, the scope block shows the frequency response
% of the system (magnitude and phase) as a function of simulation time.
%
% * In the magnitude plot, the maximum value at $t \approx 3$ indicates
% the response of the Plant block to a resonant frequency. The peak
% amplitude is approximately 1.7.
% * In the phase plot, the angle changes from 0 to $-\pi$ radians.
% Each complex pole of the system adds $-\pi/2$ radians to the phase
% angle.

sim(mdl);
open_system([mdl '/Measured Frequency Response']);

%% 
% To determine the measured resonant frequency, plot the measured magnitude
% and phase against the plant input frequency. During simulation, the model
% saves these values in a signal logging object |logsout| in the MATLAB
% workspace. You can access the logged values by using the
% |<docid:simulink_ref#buqzwj5-1 get>| method. For example, to create the
% Bode plot for the measured frequency response of the system and draw a
% red cursor at 150Hz, enter:
%
%   figure;
%   subplot(211);
%   plot(logsout.get('frequency').Values.Data, logsout.get('magnitude').Values.Data);
%   line([150 150],[0 2],'Color','red','LineStyle','--');
%   grid on;
%   title('Magnitude (Output Amplitude/Input Amplitude)');
% 
%   subplot(212);
%   plot(logsout.get('frequency').Values.Data, logsout.get('phase').Values.Data*180/pi)
%   grid on;
%   title('Phase (Degrees)');
%   close_system(mdl, 0);
%
% The Bode plot shows that the measured resonant frequency is approximately
% 150Hz, matching the value predicted by the plant dynamics.

figure;
subplot(211);
plot(logsout.get('frequency').Values.Data, logsout.get('magnitude').Values.Data);
line([150 150],[0 2],'Color','red','LineStyle','--');
grid on;
title('Magnitude (Output Amplitude/Input Amplitude)');

subplot(212);
plot(logsout.get('frequency').Values.Data, logsout.get('phase').Values.Data*180/pi)
grid on;
title('Phase (Degrees)');