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.