SDOF FRF (FFT) Magnitude discrepancy

This is my first question.

I've been having trouble with the following segment of code. For some reason, the FRF in the frequency domain is not the same magnitude as the FRF calculated by taking the FFT of an equivalent system impulse in the time domain. The phase seems equivalent, just not the magnitude.

%%%%%%%%%%%%%%%%%%SCRIPT BEGINS%%%%%%%%%%%%%%%%%

clc;

close all;

clear;

m = 6.1329; % Define mass

zeta = 0.0208; % Define damping ratio

k = 5.6367e+10; % Define stiffness

freq = [0:0.2:50000]; % Construct frequency values array

nfreq = sqrt(k/m); % Natural frequency calculation

dnfreq = sqrt(k/m)*sqrt(1-zeta^2); % Damped natural frequency calculation

%%%Constructing time array

atdl = 2*length(freq)-2; % Calculating number of elements in time array

smpprd = 1/(2*max(freq)); % Calculating data sample period

time = [0:smpprd:(atdl-1)*smpprd];% Filling time array

%%%Calculate SDOF displacement response [X(s)/F(s)] in the frequency domain.

RS = 1./(nfreq^2-(2*pi*freq).^2+j*2*zeta*(2*pi*freq)*nfreq)*(1/m);

%%%Plot the response

figure;

plot(freq,abs(RS));

hold on;

%%%Calculating the theoretical impulse response in the time domain.

xt = 1./(m*dnfreq).*exp(-zeta*nfreq*time).*sin(dnfreq*time);

%%%Getting the time domain impulse response into the frequency domain for comparison with the original frequency domain SDOF calculation

TEMP = 2*fft(xt)/atdl;

plot(freq, abs(TEMP(1:length(RS))),'r');

%%%Plot complex frequency domain impulse response for both methods 

figure;

plot3([1:atdl],real([RS,fliplr(conj(RS(2:end-1)))]),imag([RS,fliplr(conj(RS(2:end-1)))]));

hold on;

plot3([1:length(TEMP)],real(TEMP),imag(TEMP),'r');

% Calculate ratio between two methods.

n = max(abs(RS))/(max(abs(TEMP(1:length(RS)))*2/atdl))

%%%%%%%%%%%%%%%%%%%SCRIPT  ENDS%%%%%%%%%%%%%%%%%%%

Any help would be appreciated. I keep running into dead ends. I was able to figure out that the ratio between the two methods changes if you vary zeta(damping ratio).

Thanks,

Matt

Answer 1

[Dr. Seis]

I ended up copying your code and running with some modifications. Looks like it will work a little better if you make the suggestions I made above:

1. Change to TEMP = fft(xt)*smpprd;

2. Remove "2/atdl" from your bottom line (where you determine 'n')

Although, there still does appear to be some dependence on "zeta" on getting a perfect fit (i.e., the bigger zeta gets, the bigger your residuals get). Are you sure the functions you have for the time domain and frequency domain are defined properly?

Answer 2

[Matt Szelistowski]

Thank you so much for you answer.

However, I seem to be missing something. When converting a signal from the time domain to the frequency domain using FFT(), aren't the resulting magnitudes in the frequency domain supposed to be divided by the length of the original time domain and multiplied by two in order to be accurate? What is the significance of multiplying by the time domain data sample period?

[Dr. Seis]

There are two ways to explain... I will take a shot that the first way will work. Think about the forward Fourier transform for a continuous function - the function is integrated with respect to time (i.e., G(f) = integral( g(t) exp(-1i*2*pi*f*t) dt). For discrete data, the integral turns into a summation of discretized areas (length*width) under your curve defined by the amplitude length or height ( g(t) * exp(-1i*2*pi*f*t ) and the sampled width ( dt ). When Matlab does the FFT, it basically assumes dt is 1. However, this is not always the case. Therefore, the result of the FFT from Matlab needs to be multiplied by your ( dt ) in order correct the area under the curve.

When Matlab does the IFFT, it assumes that df = 1 / (N * dt), where dt is equal to 1.

Answer 3

[Matt Szelistowski]

By defined properly, do you mean are they the same SDOF system? I've mathematically derived them from each other and they are accurate as far as I know.

Answer 4

[Matt Szelistowski]

The method I used for my previous problem is illustrated by the following code. Why would calculation of the FFT() algorithm be different for the two instances?

%%%%%%%%%%%%%%%%%%%SCRIPT BEGINS%%%%%%%%%%%%%%%%%%%

clc;

close all;

clear;

smpprd = pi/20;

x = [0:smpprd:4*pi]; % time array defined

F = cos(x); % cosine function difined of magnitude 1

figure;

plot(x,F);

G = abs(fft(F)); % calculating frequency data for the signal

figure;

plot(2*G(1:end/2)/length(x)); % frequency domain data

hold on;

plot(smpprd*G(1:end/2),'r'); % using alternative method - multiplying by smpprd

% The correct output for the fft signal would have a max peak with a magnitude of 1

%%%%%%%%%%%%%%%%%%%SCRIPT ENDS%%%%%%%%%%%%%%%%%%%

[Dr. Seis]

Are you sure the correct output for the FFT signal should have a max peak with magnitude 1? In your example, the frequency of your cosine wave is 1/(2*pi). So, if in order to take the Fourier transform of a continuous function equal to cos(2*pi*ff*t), where 'ff' is the fundamental frequency of your cosine wave, on the interval from 0 to 4*pi the you would do:

ff = 1/(2*pi);

G_ff = quad(@(t)cos(2*pi*ff*t).*exp(-1i*2*pi*ff*t),0,4*pi)

which equals 6.2832 + 5.5511e-17i

When I look at your frequency spectrum plot, the peak ABS value amplitude is 6.395 for "my method" for the discrete cosine function. In fact, the smaller you make "smpprd" the closer the peak value approaches the true value obtained above.