Linear convolution using fft for system output

Question

Here is a mass-spring-damper system with an impulse response, h and an arbitrary forcing function, f (cos(t) in this case). I am trying to use Matlab's FFT function in order to perform convolution in the frequency domain. I am expecting for the output (ifft(conv)) to be the solution to the mass-spring-damper system with the specified forcing, however my plot looks completely wrong! So, i must be implementing something wrong. Please help me find my errors in my code below! Thanks

clear
%system parameters
m=4;               
k=256;                 
c=1;

wn=sqrt(k/m);
z=c/2/sqrt(m*k);
wd=wn*sqrt(1-z^2);
w=sqrt(4*k*m-c^2)/(2*m);

x0=0;   %initial conditions
v0=0;
%%%%%%%%%%%%%%%%%%%%%%%%%
t=0:.01:2*pi  ;%time vector

f=[cos(t),zeros(1,length(t)-1)];   %force f
F=fft(f);

h=[1/m/wd*exp(-z*wn*t).*sin(wd*t),zeros(1,length(t)-1)];  %impulse response
H=fft(h);

conv=H.*F;   %convolution is multiplication in freq domain

plot(0:.01:4*pi,ifft(conv))

To see what is expected run this code. Enter in cos(t); 4; 1; 256 for the inputs. You can see that it reaches a steady state at an amplitude much different than the plot generated from the above FFT code.

%%%FOR UNDERDAMPED SYSTEMS

func=input('enter function of t--->   ','s');
m=input('mass    ');
c=input('c    ');
k=input('k    ');

z=c/2/sqrt(k*m);
wn=sqrt(k/m);
wd=wn*sqrt(1-z^2);

dt=.001;
tMax=20;

x0=0;
v0=0;
t0=0;

t=0:dt:tMax;

X(:,1)=[x0;v0;t0];

functionForce=inline(func);

tic
for i=1:length(t)-1
    a=([0, 1, 0; -wn^2, -2*z*wn, 0; 0,0,0]*[X(1,i);X(2,i);X(3,i)]+[0;functionForce(X(3,i));0]);
    Xtemp=X(:,i)+[0;0;dt/2] + a*dt/2;
    b=([0, 1, 0; -wn^2, -2*z*wn, 0; 0,0,0]*[Xtemp(1);Xtemp(2);Xtemp(3)]+[0;functionForce(X(3,i));0]);
    Xtemp=Xtemp+[0;0;dt/2] + b*dt/2;
    c=([0, 1, 0; -wn^2, -2*z*wn, 0; 0,0,0]*[Xtemp(1);Xtemp(2);Xtemp(3)]+[0;functionForce(X(3,i));0]);
    Xtemp=Xtemp + [0;0;dt] +c*dt;
    d=([0, 1, 0; -wn^2, -2*z*wn, 0; 0,0,0]*[Xtemp(1);Xtemp(2);Xtemp(3)]+[0;functionForce(X(3,i));0]);
    X(:,i+1)=X(:,i)+(a+2*b+2*c+d)*dt/6+[0;0;dt];

end
toc
figure(1)
axis([0 tMax min(X(1,:)) max(X(1,:))])
plot(t,X(1,:))

What sort of plot do you expect? I get a small-amplitude/high-frequency sinusoid offset by high-amplitude/low-frequency sinusoid, which doesn't seem unreasonable for the forcing function you provided. But it's also been a while since I've studied vibrations. — eigenchris
Can you show us what you expect the output to look like? I also ran the code too and I seem to be getting something sensible like eigenchris. — rayryeng
I just added a program that will show what is expected ( i am unable to post pictures). This is a numerical solution of the system. Notice that it is WWAAAAAAAYYY too slow. The slow speed is why i am trying to use FFT method instead. Also, I'm not sure whether or not to expect to see the initial transient solution in the FFT method because of the assumption of periodicity of the input. @eigenchris — Mike James Johnson

SleuthEye SleuthEye · Accepted Answer · 2015-04-17T03:10:22

The initial transient will appear in the FFT method, so you will need to increase the time span (eg t=0:0.01:15*pi) to ensure the result includes the steady state. Incidentally since you truncate h after the same duration, increasing the time span also makes your impulse response h a better approximation to the actual infinite length impulse response.

So, updating your code to:

T=15*pi;
N=floor(T/.01);
t=[0:N-1]*0.01;  ;%time vector

...

plot([0:2*N-2]*0.01, real(ifft(conv)));
axis([0 T]); % only show the duration where the driving force is active

would correspondingly show the following response:

enter image description here

which shows the initial transient, and eventually the steady state. You may notice that the plot is similar to your alternate implementation up to a scaling factor.

This difference in scaling comes from two factors. The first one is simply due to the fact that the convolution in your FFT based implementation computes a sum which isn't weight by the time step (as compare with the dt scaling used in your second implementation). The second one comes from the fact that the second implementation does not account for the mass m for the effect of the driving force.

After accounting for those two factors, you should get the following responses:

enter image description here

where the blue curve represents the FFT based implementation (same as the above curve but scaled down by dt=0.01), and the red curve represents your alternate implementation (with the functionForce divided by m).

Linear convolution using fft for system output

1 Answers