%% Estimate the Transfer Function in MATLAB
% To estimate the transfer function of a system in MATLAB(TM), use the 
% |dsp.TransferFunctionEstimator| System object(TM). The object implements 
% the Welch's average modified periodogram method and uses the measured 
% input and output data for estimation. 

%% Initialize the System
% The system is a cascade of two filter stages: dsp.LowpassFilter and a 
% parallel connection of dsp.AllpassFilter and dsp.AllpoleFilter.
allpole = dsp.AllpoleFilter;
allpass = dsp.AllpassFilter;
lpfilter = dsp.LowpassFilter;

%% Specify Signal Source
% The input to the system is a sine wave with a frequency of 100 Hz. The 
% sampling frequency is 44.1 kHz.
sine = dsp.SineWave('Frequency',100,'SampleRate',44100,...
    'SamplesPerFrame',1024);

%% Create Transfer Function Estimator
% To estimate the transfer function of the system, create the 
% |dsp.TransferFunctionEstimator| System object. 
tfe  = dsp.TransferFunctionEstimator('FrequencyRange','onesided',...
    'OutputCoherence', true);

%% Create Array Plot
% Initialize two |dsp.ArrayPlot| objects: one to display the magnitude 
% response of the system and the other to display the coherence estimate 
% between the input and the output.
tfeplotter = dsp.ArrayPlot('PlotType','Line',...
    'XLabel','Frequency (Hz)',...
    'YLabel','Magnitude Response (dB)',...
    'YLimits',[-120 20],...
    'XOffset',0,...
    'XLabel','Frequency (Hz)',...
    'Title','System Transfer Function',...
    'SampleIncrement',44100/1024);
coherenceplotter = dsp.ArrayPlot('PlotType','Line',...
    'YLimits',[0 1.2],...
    'YLabel','Coherence',...
    'XOffset',0,...
    'XLabel','Frequency (Hz)',...
    'Title','Coherence Estimate',...
    'SampleIncrement',44100/1024);

%%
% By default, the x-axis of the array plot is in samples. To convert this 
% axis into frequency, set the 'SampleIncrement' property of the 
% |dsp.ArrayPlot| object to Fs/1024. In this example, this value is 
% 44100/1024, or 43.0664. For a two-sided spectrum, the |XOffset| property 
% of the |dsp.ArrayPlot| object must be [-Fs/2]. The frequency varies in 
% the range [-Fs/2 Fs/2]. In this example, the array plot shows a one-sided 
% spectrum. Hence, set the |XOffset| to 0. The frequency varies in the range 
% [0 Fs/2].

%% Estimate the Transfer Function
% The transfer function estimator accepts two signals: input to the two-stage 
% filter and output of the two-stage filter. The input to the filter is a 
% sine wave containing additive white Gaussian noise. The noise has a mean 
% of zero and a standard deviation of 0.1. The estimator estimates the 
% transfer function of the two-stage filter. The output of the estimator is 
% the frequency response of the filter, which is complex. To extract the 
% magnitude portion of this complex estimate, use the abs function. To convert 
% the result into dB, apply a conversion factor of 20*log10(magnitude).
for Iter = 1:1000
    input = sine() + .1*randn(1024,1);
    lpfout = lpfilter(input);
    allpoleout = allpole(lpfout);
    allpassout = allpass(lpfout);
    output = allpoleout + allpassout;
    [tfeoutput,outputcoh] = tfe(input,output);
    tfeplotter(20*log10(abs(tfeoutput)));
    coherenceplotter(outputcoh);
end
%%
% 
% <<../arrayplot_transferfunction.png>>
% 

%% 
% The first plot shows the magnitude response of the system. The second 
% plot shows the coherence estimate between the input and output of the 
% system. Coherence in the plot varies in the range [0 1] as expected.

%% Magnitude Response of the Filter Using fvtool
% The filter is a cascade of two filter stages - dsp.LowpassFilter and a 
% parallel connection of dsp.AllpassFilter and dsp.AllpoleFilter. All the 
% filter objects are used in their default state. Using the filter 
% coefficients, derive the system transfer function and plot the frequency 
% response using freqz. Below are the coefficients in the [Num] [Den] format:
% 
% * All pole filter - [1 0] [1 0.1]
% * All pass filter - [0.5 -1/sqrt(2) 1] [1 -1/sqrt(2) 0.5]
% * Lowpass filter - Determine the coefficients using the following commands:
lpf = dsp.LowpassFilter;
Coefficients = coeffs(lpf);

%%
% Coefficients.Numerator gives the coefficients in an array format. The 
% mathematical derivation of the overall system transfer function is not 
% shown here. Once you derive the transfer function, run fvtool and you can 
% see the frequency response below:
%
% <<../freqz_responseTwo.png>>
%

%%
% The magnitude response that fvtool shows matches the magnitude response 
% that the |dsp.TransferFunctionEstimator| object estimates.

% Copyright 2018 The MathWorks, Inc.