使用fft进行系统输出的线性卷积 [英] Linear convolution using fft for system output

问题描述

这是一个质量弹簧-阻尼器系统，具有脉冲响应h和任意强制功能f(在这种情况下为cos(t)).我试图使用Matlab的FFT函数以便在频域中执行卷积.我期望输出(ifft(conv))是具有指定强制的质量弹簧阻尼器系统的解决方案，但是我的绘图看起来完全错误！因此，我必须执行一些错误的操作.请帮助我在下面的代码中找到错误！谢谢

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

要查看预期结果，请运行此代码.输入cos(t)； 4; 1;输入为256.您可以看到它达到了稳态，其幅度与上述FFT代码生成的图大不相同.

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

推荐答案

初始瞬变将出现在FFT方法中，因此您需要增加时间跨度(例如t=0:0.01:15*pi)以确保结果包括稳态.顺便说一句，由于您在相同的持续时间后截断了h，因此增加时间跨度也会使您的脉冲响应h更好地逼近实际的无限长脉冲响应.

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.

缩放比例上的差异来自两个因素.第一个仅仅是因为这样一个事实，即基于FFT的实现中的卷积计算出的总和不是按时间步长加权的(与第二个实现中使用的dt缩放比例相比).第二个原因是这样的事实，第二个实现没有考虑驱动力作用的质量m.

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:

其中蓝色曲线代表基于FFT的实现(与上述曲线相同，但按dt=0.01缩小)，红色曲线代表您的替代实现(functionForce除以m).

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

这篇关于使用fft进行系统输出的线性卷积的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！