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
– eigenchris, Commented Apr 16, 2015 at 21:46
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
– rayryeng, Commented Apr 16, 2015 at 21:51
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
– Mike James Johnson, Commented Apr 17, 2015 at 0:24

SleuthEye · Accepted Answer · 2015-04-17 03:10:22Z

2

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

answered Apr 17, 2015 at 3:10

SleuthEye

14.7k2 gold badges37 silver badges67 bronze badges

Sign up to request clarification or add additional context in comments.

5 Comments

Mike James Johnson Over a year ago

ah YES! Thank you for pointing out that mistake!! I was wondering why my amplitude seemed larger than expected in both cases. One last question: How would i alter the code if i was ONLY interested in the steady state response and did not want to see the transient part? Would i need to get rid of the zero padding?

SleuthEye Over a year ago

No, you need the zero padding (otherwise you'd get a circular convolution). If you want the steady-state part, get rid of the initial samples of the convolution results (e.g. convresult = real(ifft(conv)); result = convresult(delay+1:end);, to get rid of the first delay samples).

Mike James Johnson Over a year ago

That makes sense; however, i was hoping there is a way to not have to calcualte which part of the signal to throw out, but rather have that taken care of through the FFT calculations. Maybe there can be an advantage to using circular convolution in order to ignore the transient part? If there is a force with a defined period (cos(t) ; 2*pi ), and only the steady state is of interest, wouldn't circular convolution be preferred over linear?

SleuthEye Over a year ago

Main problem with that is that the impulse response is not periodic, so circular convolution would't make sense in that case. I'd look at the possibility to reformulate the equations with periodic boundary conditions to avoid using the impulse response.

Mike James Johnson Over a year ago

Ok, your responses are much appreciated! I will work on that. If i run into trouble I may post in another question.

Collectives™ on Stack Overflow

Linear convolution using fft for system output

1 Answer 1

5 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

5 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related