SDOF FRF (FFT) Magnitude discrepancy

18 views (last 30 days)
Matt Szelistowski on 6 Jan 2012
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

Dr. Seis on 6 Jan 2012
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?

Matt Szelistowski on 6 Jan 2012
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 on 7 Jan 2012
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.

Matt Szelistowski on 6 Jan 2012
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.

Matt Szelistowski on 6 Jan 2012
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 on 7 Jan 2012
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);