www.gusucode.com > dsp 案例源码程序 matlab代码 > dsp/EstimateTheTransferFunctionInMATLABExample.m

    %% 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 $\frac{F_{s}}{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

%%
% 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:
%
% <<../freqz_response.png>>

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