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,:))

Answer 1

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:

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:

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).