Why am I getting NAN values for magnitudes when I use tfestimate?

Hi all,

I am simulating a high frequency (6250 Hz) mass, spring and damper system's response to an impulse input.

My objective is to estimate how much noise can be on the measurement output to be able to determine the TF of the system with 1% accuracy in magnitude and 0.1% accuracy in frequency.

For this objective I am using the following system parameters:

% system parameters in the example
m = 9; %kg
d = 0.1; % damper value
res_freq = (6250)*2*pi; %6250 Hz I want to have a resonance peak.
k = m*(res_freq)^2; %stiffness value
s =tf('s');
G = 1/(m*s^2 + d*s +k);
Gd = c2d(G,5e-5);
opts = bodeoptions('cstprefs');
opts.PhaseVisible = 'off';
opts.FreqUnits = 'Hz';
figure;
bode(Gd,opts);
hold on
bode(G,opts);
grid on;

with the lines on plotting the bode of the continuous and discrete systems, I am able to confirm that the example is a valid one.

The system is simulated as follows:

Fs = 20000; %sampling rate is 20000 Hz
Ts = 1/Fs;
Ttot = [0: Ts: (2-Ts)]';

% create an impulse input for the length of 2 seconds.
L = length(Ttot);
dirac_input = 200 * ones(1,1); % dirac impulse with 200N force magnitude
input = [dirac_input; zeros(39999,1)];

x0 = [0 0];
oupt_NoNoise = lsim(Gd,input,Ttot,x0); %output from impulse without noise

%Adding noise to the simulation value: Tune this expression to understand noise sensitvity...
oupt_Noise = oupt_NoNoise+ (max(oupt_NoNoise)/100)*randn(length(oupt_NoNoise),1);  % adding 1% noise.

With the output with and without noise simulated I performed the FFT of the two signals and I do notice the peak on the resonance frequency of the system as the most dominant freq. The signal with the noise had small amplitudes of other frequencies as well. So the output data is also as per my expectation.

Now, I am trying to see if I can recreate the system TF based on input and output data. So I do the following:

%% TF estimate with data that is noise free as well as data that is noisy.
nfft = 16384; 
noverlap = nfft/2;
window = hann(nfft);
[txy,f] = tfestimate(input,oupt_NoNoise,window,noverlap,[],Fs); %getting tf estimate on data with no noise.
[txy_n,f_n] = tfestimate(input,oupt_Noise,window,noverlap,[],Fs); %getting tf estimate on data with noise.

figure;
plot(f, 20*log10(abs(txy)),'r');  
hold on;
plot(f_n, 20*log10(abs(txy_n)),'b'); 
grid on;

and when I do this, I get that the values in variables txy and txy_n to be NAN. What is it that I am doing wrong in terms of using the tfestimate, can anybody please advice?

[Wayne King]

Hi Harryboy, Can you please provide the input vector? I don't see input defined here.

[Harryboy]

Hi Wayne King, thanks for the comment. I have updated the input in my queston.

to make it easy for you, I am pasting it below as well:

% create an impulse input for the length of 2 seconds.
L = length(Ttot);
dirac_input = 200 * ones(1,1); % dirac impulse with 200N force magnitude
input = [dirac_input; zeros(39999,1)];

FYI:

You may indeed notice, that I create an impulse input vector. For simulating the output I choose the inbuilt matlab command 'impulse' rather than using lsim. The size of the input and the output vector are the same since this is a SISO system.

When I did try to use the lsim command to simulate the output I get a warning that the system is undersampled, which is an incorrect statement in my opinion.

Since the sizes remain same, I would expect that using 'impulse' is still ok to get the o/p. And then using the output and the input vector that was created can be used to estimate the TF. Do you have any comments on this implementation?

[Harryboy]

Hi,

when I change the input from:

dirac_input = 200 * ones(1,1); % ideal impulse.
input = [dirac_input; zeros(39999,1)]; 

to simulating a step:

dirac_input = 200 * ones(1,1); % ideal impulse.
input = [dirac_input; 200*ones(39999,1)]; % I have made it into a step so that I have output without NAN in the tfestimate output

then I am able to get non NAN output from the tfestimate command. With the impulse then it is clearly diving something with zero. But I imagined that all the division is happening at frequency domain and an ideal impulse has non zero magnitude at all frequencies.

Anyone able to explain this observation?