%% Group Delay Estimation
%%
% *Note*: This example runs only in R2016b or later. If you are using an
% earlier release, replace each call to the function with the equivalent
% |step| syntax. For example, myObject(x) becomes step(myObject,x).

%%
% Estimate the group delay of a linear phase FIR filter using a 
% |dsp.TransferFunctionEstimator| object followed by |dsp.PhaseExtractor| 
% and |dsp.Differentiator| objects. The group delay of a linear phase FIR
% filter is given by $GD = -(d\theta(\omega)/d\omega) = - \frac{N}{2}$, 
% where $\theta(\omega)$ is the phase information of the filter, $\omega$ is the 
% frequency vector, and _N_ is the order of the filter.

%% Set Up the Objects
% Create a linear phase FIR lowpass filter. Set the order to 200, the 
% passband frequency to 255 Hz, the passband ripple to 0.1 dB, and the 
% stopband attenuation to 80 dB. Specify a sample rate of 512 Hz.
Fs = 512;
LPF = dsp.LowpassFilter('SampleRate',Fs,'PassbandFrequency',255,...
    'DesignForMinimumOrder',false,'FilterOrder',200);
%%
% To estimate the transfer function of the lowpass filter, create a 
% transfer function estimator. Specify the window to be |Hann|.
% Set the FFT length to 1024 and the number of spectral averages to
% 200. 
TFE = dsp.TransferFunctionEstimator('FrequencyRange','twosided',...
    'SpectralAverages',200,'FFTLengthSource','Property',...
    'FFTLength',1024);
%%
% To extract the unwrapped phase from the frequency response of the filter,
% create a phase extractor.
PE = dsp.PhaseExtractor;
%%
% To differentiate the phase $\theta$, create a differentiator filter.  
% This value is used in computing the group delay.
DF = dsp.Differentiator;

%%
% To smoothen the input, create a variable bandwidth FIR filter.
Gain1 = 512/pi;
Gain2 = -1;
VBFilter = dsp.VariableBandwidthFIRFilter('CutoffFrequency',10,...
    'SampleRate',Fs);

%%
% To view the group delay of the filter, create an array plot object.
AP = dsp.ArrayPlot('PlotType','Line','YLimits',[-500 400],...
    'YLabel','Amplitude','XLabel','Number of samples');

%% Run the Algorithm
% The |for|-loop is the streaming loop that estimates the group delay of the 
% filter. In the loop, the algorithm filters the input signal, estimates 
% the transfer function of the filter, and differentiates the phase of the 
% filter to compute the group delay.
Niter = 1000; % Number of iterations
for k = 1:Niter
        x = randn(512,1);  % Input signal = white Gaussian noise
        y = LPF(x);   % Filter noise with Lowpass FIR filter
        H = TFE(x,y); % Compute transfer function estimate
        Phase = PE(H); % Extract the Unwrapped phase
        phaseaftergain1 = Gain1*Phase;
        DiffOut = DF(phaseaftergain1); % Differentiate the phase
        phaseaftergain2 = Gain2 * DiffOut; 
        VBFOut = VBFilter(phaseaftergain2); % Smooth the group delay
        AP(VBFOut); % Display the group delay
end
%%
% 
% <<../arrayplot_groupdelayestimation.png>>
% 

%%
% As you can see, the group delay of the lowpass filter is 100.

%% 
% Copyright 2012 The MathWorks, Inc.