www.gusucode.com > ident_featured 案例代码 matlab源码程序 > ident_featured/iddemofr.m
%% Frequency Domain Identification: Estimating Models Using Frequency Domain Data % This example shows how to estimate models using frequency domain data. % The estimation and validation of models using frequency domain data work % the same way as they do with time domain data. This provides a great % amount of flexibility in estimation and analysis of models using time and % frequency domain as well as spectral (FRF) data. You may simultaneously % estimate models using data in both domains, compare and combine these % models. A model estimated using time domain data may be validated using % spectral data or vice-versa. % % Frequency domain data can not be used for estimation or validation of % nonlinear models. % Copyright 1986-2015 The MathWorks, Inc. %% Introduction % Frequency domain experimental data are common in many applications. It % could be that the data was collected as frequency response data % (frequency functions: FRF) from the process using a frequency analyzer. % It could also be that it is more practical to work with the input's and % output's Fourier transforms (FFT of time-domain data), for example to % handle periodic or band-limited data. (A band-limited continuous time % signal has no frequency components above the Nyquist frequency). In % System Identification Toolbox, frequency domain I/O data are represented % the same way as time-domain data, i.e., using |iddata| objects. The % 'Domain' property of the object must be set to 'Frequency'. Frequency % response data are represented as complex vectors or as magnitude/phase % vectors as a function of frequency. IDFRD objects in the toolbox are used % to encapsulate FRFs, where a user specifies the complex response data and % a frequency vector. Such IDDATA or IDFRD objects (and also FRD objects of % Control System Toolbox) may be used seamlessly with any estimation % routine (such as |procest|, |tfest| etc). %% Inspecting Frequency Domain Data % Let us begin by loading some frequency domain data: load demofr %% % This MAT-file contains frequency response data at frequencies |W|, with % the amplitude response |AMP| and the phase response |PHA|. Let us first % have a look at the data: subplot(211), loglog(W,AMP),title('Amplitude Response') subplot(212), semilogx(W,PHA),title('Phase Response') %% % This experimental data will now be stored as an IDFRD object. First % transform amplitude and phase to a complex valued response: zfr = AMP.*exp(1i*PHA*pi/180); Ts = 0.1; gfr = idfrd(zfr,W,Ts); %% % |Ts| is the sample time of the underlying data. If the data % corresponds to continuous time, for example since the input has been % band-limited, use Ts = 0. % % Note: If you have the Control System Toolbox(TM), you could use an FRD % object instead of the IDFRD object. IDFRD has options for more % information, like disturbance spectra and uncertainty measures which are % not available in FRD objects. % % The IDFRD object |gfr| now contains the data, and it can be plotted and % analyzed in different ways. To view the data, we may use |plot| or % |bode|: clf bode(gfr), legend('gfr') %% Estimating Models Using Frequency Response (FRF) Data % To estimate models, you can now use |gfr| as a data set with all the % commands of the toolbox in a transparent fashion. The only restriction % is that noise models cannot be built. This means that for polynomial % models only OE (output-error models) apply, and for state-space models, % you have to fix |K = 0|. m1 = oe(gfr,[2 2 1]) % Discrete-time Output error (transfer function) model ms = ssest(gfr) % Continuous-time state-space model with default choice of order mproc = procest(gfr,'P2UDZ') % 2nd-order, continuous-time model with underdamped poles compare(gfr,m1,ms,mproc) L = findobj(gcf,'type','legend'); L.Location = 'southwest'; % move legend to non-overlapping location %% % As shown above a variety of linear model types may be estimated in both % continuous and discrete time domains, using spectral data. These models % may be validated using, time-domain data. The time-domain I/O data set % |ztime|, for example, is collected from the same system, and can be used % for validation of |m1|, |ms| and |mproc|: compare(ztime,m1,ms,mproc) %validation in a different domain %% % We may also look at the residuals to affirm the quality of the model % using the validation data |ztime|. As observed, the residuals are almost % white: resid(ztime,mproc) % Residuals plot %% Condensing Data Using SPAFDR % An important reason to work with frequency response data is that it is % easy to condense the information with little loss. The command SPAFDR % allows you to compute smoothed response data over limited frequencies, % for example with logarithmic spacing. Here is an example where the |gfr| % data is condensed to 100 logarithmically spaced frequency values. With % a similar technique, also the original time domain data can be % condensed: sgfr = spafdr(gfr) % spectral estimation with frequency-dependent resolution sz = spafdr(ztime); % spectral estimation using time-domain data clf bode(gfr,sgfr,sz) axis([pi/100 10*pi, -272 105]) legend('gfr (raw data)','sgfr','sz','location','southwest') %% % The Bode plots show that the information in the smoothed data has been % taken well care of. Now, these data records with 100 points can very well % be used for model estimation. For example: msm = oe(sgfr,[2 2 1]); compare(ztime,msm,m1) % msm has the same accuracy as M1 (based on 1000 points) %% Estimation Using Frequency-Domain I/O Data % It may be that the measurements are available as Fourier transforms of % inputs and output. Such frequency domain data from the system are given % as the signals Y and U. In loglog plots they look like Wfd = (0:500)'*10*pi/500; subplot(211),loglog(Wfd,abs(Y)),title('The amplitude of the output') subplot(212),loglog(Wfd,abs(U)),title('The amplitude of the input') %% % The frequency response data is essentially the ratio between |Y| and |U|. % To collect the frequency domain data as an IDDATA object, do as follows: ZFD = iddata(Y, U, 'Ts', 0.1, 'Frequency', Wfd) %% % Now, again the frequency domain data set |ZFD| can be used as data in all % estimation routines, just as time domain data and frequency response % data: mf = ssest(ZFD) % SSEST picks best order in 1:10 range when called this way mfr = ssregest(ZFD) % an alternative regularized reduction based state-space estimator clf compare(ztime,mf,mfr,m1) %% Transformations Between Data Representations (Time - Frequency) % Time and frequency domain input-output data sets can be transformed to % either domain by using FFT and IFFT. These commands are adapted to % IDDATA objects: dataf = fft(ztime) datat = ifft(dataf) %% % Time and frequency domain input-output data can be transformed to % frequency response data by SPAFDR, SPA and ETFE: g1 = spafdr(ztime) g2 = spafdr(ZFD); clf; bode(g1,g2) %% % Frequency response data can also be transformed to more smoothed data % (less resolution and less data) by SPAFDR and SPA; g3 = spafdr(gfr); %% % Frequency response data can be transformed to frequency domain % input-output signals by the command IDDATA: gfd = iddata(g3) %% Using Continuous-time Frequency-domain Data to Estimate Continuous-time Models % Time domain data can naturally only be stored and dealt with as % discrete-time, sampled data. Frequency domain data have the advantage % that continuous time data can be represented correctly. Suppose that the % underlying continuous time signals have no frequency information above % the Nyquist frequency, e.g. because they are sampled fast, or the input % has no frequency component above the Nyquist frequency. Then the Discrete % Fourier transforms (DFT) of the data also are the Fourier transforms of % the continuous time signals, at the chosen frequencies. They can % therefore be used to directly fit continuous time models. In fact, this % is the correct way of using band-limited data for model fit. % % This will be illustrated by the following example. %% % Consider the continuous time system: % % $$ G(s) = \frac{1}{s^2+s+1} $$ % m0 = idpoly(1,1,1,1,[1 1 1],'ts',0) %% % Choose an input with low frequency contents that is fast sampled: rng(235,'twister'); u = idinput(500,'sine',[0 0.2]); u = iddata([],u,0.1,'intersamp','bl'); %% % 0.1 is the sample time, and |'bl'| indicates that the input is % band-limited, i.e. in continuous time it consists of sinusoids % with frequencies below half the sampling frequency. Correct simulation of % such a system should be done in the frequency domain: uf = fft(u); uf.ts = 0; % Denoting that the data is continuous time yf = sim(m0,uf); % % Add some noise to the data: yf.y = yf.y + 0.05*(randn(size(yf.y))+1i*randn(size(yf.y))); dataf = [yf uf] % This is now a continuous time frequency domain data set. %% % Look at the data: plot(dataf) axis([0 10 0 2.5]) %% % Using |dataf| for estimation will by default give continuous time models: % State-space: m4 = ssest(dataf,2); %Second order continuous-time model %% % For a polynomial model with |nb = 2| numerator coefficient and |nf = 2| % estimated denominator coefficients use: nb = 2; nf = 2; m5 = oe(dataf,[nb nf]) %% % Compare step responses with uncertainty of the true system |m0| and the % models |m4| and |m5|. The confidence intervals are shown with patches in % the plot. clf h = stepplot(m0,m4,m5); showConfidence(h,1) legend('show','location','southeast') %% % Although it was not necessary in this case, it is generally advised to % focus the fit to a limited frequency band (low pass filter the data) % when estimating using continuous time data. The system has a bandwidth % of about 3 rad/s, and was excited by sinusoids up to 6.2 rad/s. % A reasonable frequency range to focus the fit to is then [0 7] rad/s: % m6 = ssest(dataf,2,ssestOptions('WeightingFilter',[0 7])) % state space model %% m7 = oe(dataf,[1 2],oeOptions('WeightingFilter',[0 7])) % polynomial model of Output Error structure %% opt = procestOptions('SearchMethod','lsqnonlin',... 'WeightingFilter',[0 7]); % Requires Optimization Toolbox(TM) m8 = procest(dataf,'P2UZ',opt) % process model with underdamped poles %% opt = tfestOptions('SearchMethod','lsqnonlin',... 'WeightingFilter',[0 7]); % Requires Optimization Toolbox m9 = tfest(dataf,2,opt) % transfer function with 2 poles %% h = stepplot(m0,m6,m7,m8,m9); showConfidence(h,1) legend('show') %% Conclusions % We saw how time, frequency and spectral data can seamlessly be used to % estimate a variety of linear models in both continuous and discrete time % domains. The models may be validated and compared in domains different % from the ones they were estimated in. The data formats (time, frequency % and spectrum) are interconvertible, using methods such as |fft|, |ifft|, % |spafdr| and |spa|. Furthermore, direct, continuous-time estimation is % achievable by using |tfest|, |ssest| and |procest| estimation routines. % The seamless use of data in any domain for estimation and analysis is an % important feature of System Identification Toolbox.