MATLAB 使用 FFT 找到对周期性输入力的稳态响应（质量 spring 阻尼系统）
MATLAB using FFT to find steady state response to a periodic input force (mass spring damper system)
假设我有一个质量-spring-阻尼器系统...
这是我的代码 (matlab)...
% system parameters
m=4; k=256; c=1; wn=sqrt(k/m); z=c/2/sqrt(m*k); wd=wn*sqrt(1-z^2);
% initial conditions
x0=0; v0=0;
%% time
dt=.001; tMax=2*pi; t=0:dt:tMax;
% input
F=cos(t); Fw=fft(F);
% impulse response function
h=1/m/wd*exp(-z*wn*t).*sin(wd*t); H=fft(h);
% convolution
convolution=Fw.*H; sol=ifft(convolution);
% plot
plot(t,sol)
所以我可以成功地检索一个绘图，但是我得到了奇怪的响应我还编写了一个 RK4 方法来求解微分方程组所以我知道绘图应该是什么样子，以及我从使用中得到的绘图FFT 的幅度为 2，而幅度应为 0.05。
那么，我如何使用 FFT 求解该系统的稳态响应。我想使用 FFT，因为它比数值积分方法快大约 3 个数量级。
请记住，我将我的周期性输入定义为 cos(t)，周期为 2*pi，这就是为什么我只在跨越 0 到 2*pi（1 个周期）的时间向量上使用 FFT .我还注意到，如果我将 tMax 时间更改为 2*pi 的倍数，例如 10*pi，我会得到一个看起来类似的图，但振幅是 4 而不是 2，无论哪种方式仍然不是 .05！。也许我需要乘以某种因素？
我也绘制了：plot(t,Fw) 由于强制函数是 cos(t)，我期望在 1 处看到一个峰值，但我没有看到任何峰值（也许我不应该绘制 Fw 与 t)
我知道可以使用傅立叶变换/fft 求解稳态响应，我只是漏掉了一些东西！我需要帮助和理解！！
原始结果
运行 您提供的代码并将结果与发布在 中的 RK4 代码进行比较，我们得到以下响应：
其中蓝色曲线代表基于 FFT 的实现，红色曲线代表您的备用 RK4 实现。正如你所指出的，曲线是完全不同的。
得到正确的回应
最明显的问题当然是振幅，这个问题贴出的代码振幅差异的主要来源和我在中指出的相同：
- RK4 实现执行数值积分，通过积分步骤
dt 正确缩放求和值。基于 FFT 的实现缺少这种缩放。
- 基于 FFT 的实现中使用的脉冲响应与按质量缩放的驱动力一致
m，这是 RK4 实现中缺少的一个因素。
解决这两个问题会导致响应更接近一些，但仍然不完全相同。正如您可能发现的那样，鉴于您的其他问题的已发布代码发生了变化，另一件缺少的事情是输入和脉冲响应的零填充，没有它你会得到一个 circular 卷积而不是线性卷积：
f=[cos(t),zeros(1,length(t)-1)]; %force f
h=[1/m/wd*exp(-z*wn*t).*sin(wd*t),zeros(1,length(t)-1)]; %impulse response
最后，确保卷积产生良好结果的最后一个要素是对无限长度脉冲响应使用良好的近似。多长时间足够长取决于脉冲响应的衰减率。使用您提供的参数，脉冲响应会在大约 11*pi 后减弱到其原始值的 1%。因此，将时间跨度扩展到 tMax=14*pi（在脉冲响应消失后包括一个完整的 2*pi 周期）可能就足够了。
正在获取 steady-state 响应
获得 steady-state 响应的最简单方法是丢弃初始瞬态。在这个过程中我们丢弃了整数个周期的参考驱动力（当然这需要知道驱动力的基频）：
T0 = tMax-2*pi;
delay = min(find(t>T0));
sol = sol(delay:end);
plot([0:length(sol)-1]*dt, sol, 'b');
axis([0 2*pi]);
然后得到的响应是：
蓝色曲线再次代表基于 FFT 的实现，红色曲线代表您的替代 RK4 实现。好多了！
另一种方法
计算等待瞬态响应消失的许多周期的响应并提取相应的剩余样本
到稳定状态可能看起来有点浪费，尽管由于 FFT，计算仍然相当快。
那么，让我们稍微回头看看问题域。正如你可能知道的那样，
mass-spring-damper system 由微分方程控制：
其中 f(t) 是这种情况下的驱动力。
请注意，齐次方程的通解具有以下形式：
然后关键是要认识到c>0和m>0在稳定状态下消失（t趋于无穷）的情况下的通解。
因此 steady-state 解仅依赖于 non-homogenous 方程的特定解。
对于形式为
的驱动力，可以通过 method of undetermined coefficients 找到此特定解
相应地假设解的形式为
代入微分方程得到方程：
因此，解决方案可以实现为：
EF0 = [wn*wn-w*w 2*z*wn*w; -2*z*wn*w wn*wn-w*w]\[1/m; 0];
sol = EF0(1)*cos(w*t)+EF0(2)*sin(w*t);
plot(t, sol);
你的情况 w=2*pi。
泛化
通过将驱动力表示为
傅里叶级数（假设驱动力函数满足Dirichlet conditions）：
特定解可以相应地假设为
可以用与前面的情况非常相似的方式来解决特定的解决方案。这导致以下实施：
% normalize
y = F/m;
% compute coefficients proportional to the Fourier series coefficients
Yw = fft(y);
% setup the equations to solve the particular solution of the differential equation
% by the method of undetermined coefficients
k = [0:N/2];
w = 2*pi*k/T;
A = wn*wn-w.*w;
B = 2*z*wn*w;
% solve the equation [A B;-B A][real(Xw); imag(Xw)] = [real(Yw); imag(Yw)] equation
% Note that solution can be obtained by writing [A B;-B A] as a scaling + rotation
% of a 2D vector, which we solve using complex number algebra
C = sqrt(A.*A+B.*B);
theta = acos(A./C);
Ywp = exp(j*theta)./C.*Yw([1:N/2+1]);
% build a hermitian-symmetric spectrum
Xw = [Ywp conj(fliplr(Ywp(2:end-1)))];
% bring back to time-domain (function synthesis from Fourier Series coefficients)
x = ifft(Xw);
最后的说明
我在上面的推导中特意避开了无阻尼的c=0情况。在这种情况下，振荡永远不会消失，齐次方程的通解不一定是平凡的。
在这种情况下，最终的"steady-state"可能与驱动力具有相同的周期，也可能不具有相同的周期。事实上，如果通解的周期振荡与有理数（整数比）的驱动力周期无关，则它可能根本不是周期性的。
假设我有一个质量-spring-阻尼器系统...
这是我的代码 (matlab)...
% system parameters
m=4; k=256; c=1; wn=sqrt(k/m); z=c/2/sqrt(m*k); wd=wn*sqrt(1-z^2);
% initial conditions
x0=0; v0=0;
%% time
dt=.001; tMax=2*pi; t=0:dt:tMax;
% input
F=cos(t); Fw=fft(F);
% impulse response function
h=1/m/wd*exp(-z*wn*t).*sin(wd*t); H=fft(h);
% convolution
convolution=Fw.*H; sol=ifft(convolution);
% plot
plot(t,sol)
所以我可以成功地检索一个绘图，但是我得到了奇怪的响应我还编写了一个 RK4 方法来求解微分方程组所以我知道绘图应该是什么样子，以及我从使用中得到的绘图FFT 的幅度为 2，而幅度应为 0.05。
那么，我如何使用 FFT 求解该系统的稳态响应。我想使用 FFT，因为它比数值积分方法快大约 3 个数量级。
请记住，我将我的周期性输入定义为 cos(t)，周期为 2*pi，这就是为什么我只在跨越 0 到 2*pi（1 个周期）的时间向量上使用 FFT .我还注意到，如果我将 tMax 时间更改为 2*pi 的倍数，例如 10*pi，我会得到一个看起来类似的图，但振幅是 4 而不是 2，无论哪种方式仍然不是 .05！。也许我需要乘以某种因素？
我也绘制了：plot(t,Fw) 由于强制函数是 cos(t)，我期望在 1 处看到一个峰值，但我没有看到任何峰值（也许我不应该绘制 Fw 与 t)
我知道可以使用傅立叶变换/fft 求解稳态响应，我只是漏掉了一些东西！我需要帮助和理解！！
原始结果
运行 您提供的代码并将结果与发布在
其中蓝色曲线代表基于 FFT 的实现，红色曲线代表您的备用 RK4 实现。正如你所指出的，曲线是完全不同的。
得到正确的回应
最明显的问题当然是振幅，这个问题贴出的代码振幅差异的主要来源和我在
- RK4 实现执行数值积分，通过积分步骤
dt正确缩放求和值。基于 FFT 的实现缺少这种缩放。 - 基于 FFT 的实现中使用的脉冲响应与按质量缩放的驱动力一致
m，这是 RK4 实现中缺少的一个因素。
解决这两个问题会导致响应更接近一些，但仍然不完全相同。正如您可能发现的那样，鉴于您的其他问题的已发布代码发生了变化，另一件缺少的事情是输入和脉冲响应的零填充，没有它你会得到一个 circular 卷积而不是线性卷积：
f=[cos(t),zeros(1,length(t)-1)]; %force f
h=[1/m/wd*exp(-z*wn*t).*sin(wd*t),zeros(1,length(t)-1)]; %impulse response
最后，确保卷积产生良好结果的最后一个要素是对无限长度脉冲响应使用良好的近似。多长时间足够长取决于脉冲响应的衰减率。使用您提供的参数，脉冲响应会在大约 11*pi 后减弱到其原始值的 1%。因此，将时间跨度扩展到 tMax=14*pi（在脉冲响应消失后包括一个完整的 2*pi 周期）可能就足够了。
正在获取 steady-state 响应
获得 steady-state 响应的最简单方法是丢弃初始瞬态。在这个过程中我们丢弃了整数个周期的参考驱动力（当然这需要知道驱动力的基频）：
T0 = tMax-2*pi;
delay = min(find(t>T0));
sol = sol(delay:end);
plot([0:length(sol)-1]*dt, sol, 'b');
axis([0 2*pi]);
然后得到的响应是：
蓝色曲线再次代表基于 FFT 的实现，红色曲线代表您的替代 RK4 实现。好多了！
另一种方法
计算等待瞬态响应消失的许多周期的响应并提取相应的剩余样本 到稳定状态可能看起来有点浪费，尽管由于 FFT，计算仍然相当快。
那么，让我们稍微回头看看问题域。正如你可能知道的那样， mass-spring-damper system 由微分方程控制：
其中 f(t) 是这种情况下的驱动力。
请注意，齐次方程的通解具有以下形式：
然后关键是要认识到c>0和m>0在稳定状态下消失（t趋于无穷）的情况下的通解。
因此 steady-state 解仅依赖于 non-homogenous 方程的特定解。
对于形式为
相应地假设解的形式为
代入微分方程得到方程：
因此，解决方案可以实现为：
EF0 = [wn*wn-w*w 2*z*wn*w; -2*z*wn*w wn*wn-w*w]\[1/m; 0];
sol = EF0(1)*cos(w*t)+EF0(2)*sin(w*t);
plot(t, sol);
你的情况 w=2*pi。
泛化
通过将驱动力表示为 傅里叶级数（假设驱动力函数满足Dirichlet conditions）：
特定解可以相应地假设为
可以用与前面的情况非常相似的方式来解决特定的解决方案。这导致以下实施：
% normalize
y = F/m;
% compute coefficients proportional to the Fourier series coefficients
Yw = fft(y);
% setup the equations to solve the particular solution of the differential equation
% by the method of undetermined coefficients
k = [0:N/2];
w = 2*pi*k/T;
A = wn*wn-w.*w;
B = 2*z*wn*w;
% solve the equation [A B;-B A][real(Xw); imag(Xw)] = [real(Yw); imag(Yw)] equation
% Note that solution can be obtained by writing [A B;-B A] as a scaling + rotation
% of a 2D vector, which we solve using complex number algebra
C = sqrt(A.*A+B.*B);
theta = acos(A./C);
Ywp = exp(j*theta)./C.*Yw([1:N/2+1]);
% build a hermitian-symmetric spectrum
Xw = [Ywp conj(fliplr(Ywp(2:end-1)))];
% bring back to time-domain (function synthesis from Fourier Series coefficients)
x = ifft(Xw);
最后的说明
我在上面的推导中特意避开了无阻尼的c=0情况。在这种情况下，振荡永远不会消失，齐次方程的通解不一定是平凡的。
在这种情况下，最终的"steady-state"可能与驱动力具有相同的周期，也可能不具有相同的周期。事实上，如果通解的周期振荡与有理数（整数比）的驱动力周期无关，则它可能根本不是周期性的。