系统输出使用 fft 的线性卷积

Question

这是一个质量-spring-阻尼器系统，具有脉冲响应 h 和任意力函数 f（在本例中为 cos(t)）。我正在尝试使用 Matlab 的 FFT 函数在频域中执行卷积。我期望输出 (ifft(conv)) 是具有指定强制的质量-spring-阻尼器系统的解决方案，但是我的情节看起来完全错误！所以，我一定是在实施错误的东西。请帮助我在下面的代码中找到我的错误！谢谢

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 代码生成的图大不相同的幅度达到稳定状态。

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

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

因此，将您的代码更新为：

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

会相应地显示以下响应：

显示初始瞬态，最终稳定状态。您可能会注意到该图与您的替代实施方案相似，但比例因子不同。

这种缩放差异来自两个因素。第一个仅仅是因为基于 FFT 的实现中的卷积计算的总和不是时间步长的权重（与第二个实现中使用的 dt 缩放相比）。第二个来自于这样一个事实，即第二个实现没有考虑质量 m 对驱动力的影响。

考虑到这两个因素后，您应该得到以下回复：

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

系统输出使用 fft 的线性卷积

Linear convolution using fft for system output

matlab

fft

linear-algebra

convolution