Estimate the Transfer Function in MATLAB
To estimate the transfer function of a system in MATLAB, use the dsp.TransferFunctionEstimator System object. 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 . 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
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 mathemtical 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:
The magnitude response that fvtool shows matches the magnitude response that the dsp.TransferFunctionEstimator object estimates.