How to show the power spectral density of a time-domain arx model in System Identification Toolbox?

Hey all,

After constructing my time domain arx model in System Identification Toolbox, I want to transform the predicted data from the model to frequency in order to examine its psd, here is my code:

[num] = xlsread( sprintf('%d.xlsx',1),1); RRI = num(:,3); JTI = num(:,10); SamplesPerWindow = 60; ts = 0.04; % Sampling Interval fs=1/ts;

for i = 1:5

    u = RRI(((SamplesPerWindow*(i-1))+1):SamplesPerWindow*i,1);
    y = JTI(((SamplesPerWindow*(i-1))+1):SamplesPerWindow*i,1);
    data = iddata(y,u,ts); 

    advice(data); % 

    datae = misdata(data,idpoly); 
    advice(datae);

    estimated_delay = delayest(datae) %computes the delay from the input to the output
    NN = struc(2,1:20,1:20); % according to advice the order of the model cannot exceed 20 
    datae_e = datae(1:30);
    datae_v = datae(31:60);
    V = arxstruc(datae_e, datae_v, NN);
    nn = selstruc(V,'aic'); % c= aic or mdl or zero which indicates the best fit
    m = arx(datae,nn);

    %model = arx(datae, [nn(1,1) nn(1,2) estimated_delay]); % 70% = [2 9 1] , 79% =[2 15 7] , 86% = [2 20 7]
    model2 = nlarx(datae, nn, 'linear'); %71% = [2 9 1], 80.82% =[2 15 7] , 87% = [2 20 7]
    %model3 = nlarx(datae, [nn(1,1) nn(1,2) estimated_delay], 'linear');
    model4 = oe(datae,nn)
    figure, compare(datae, m, model2, model4) 

    MODEL = fft(m);
    N=length(u);
    fi=fs/N;
    f=[0:N-1]*fi;
    plot(f,abs(X))

end

Now, I'm getting an error: ??? Undefined function or method 'fft' for input arguments of type 'idpoly'.

I know that fft works for iddata objects but how can I make it work for idpoly and idnlarx???

Is it enough to write direct ffplot(m) without making any other steps after constructing the model?? Or I should first use spa command first?? I tried using it (g = spa(m)) but I am having an error: The value of the "OutputData" property must be a double matrix or a cell array of such matrices. Error in ==> spa at 56 data = iddata(data(:,1:ny),data(:,ny+1:end));

Answer 1

[Rajiv Singh]

For frequency response of the model (H(w) = B(w)/A(w)), use FREQRESP or BODE command, as in: H = freqresp(m, w); or [Magnitude, Phase, w] = bode(m);

For the disturbance spectra (PHI(w) = L*||A(w)||^2), first extract the noise model and then call the spectra computing command:

• In releases R2011b and earlier, do mn = m(:,'noise'); DS = freqresp(mn)
• In releases R2012a and later, do: mn = noise2meas(m); DS = spectrum(mn)

Frequency response for a nonlinear model is not well-defined. You can linearize the model about a given operating point (see idnlarx/linearize) and call BODE/FREQRESP on the resulting linear model.

Or, you can construct the response of the nonlinear model over a grid of frequencies by exciting the system with one sinusoid at a time. This is more easily done in Simulink if you import the identified model using the Nonlinear ARX block. Then, use FRESTIMATE (requires Simulink Control Design).

Or, if you are able to write the governing equations for the model, you can do a (painstaking) describing function/harmonic balance analysis by hand to compute the "frequency response".