MATLAB Answers

0

SDOF FRF (FFT) Magnitude discrepancy

Asked by 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

  0 Comments

Sign in to comment.

4 Answers

Answer by 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?

  0 Comments

Sign in to comment.


Answer by 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?

  1 Comment

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.

Sign in to comment.


Answer by 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.

  0 Comments

Sign in to comment.


Answer by 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%%%%%%%%%%%%%%%%%%%

  1 Comment

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.

Sign in to comment.