FFT：当信号不是"straight"时提取振幅比
FFT: extract amplitude ratio when signal is not "straight"
我需要对一个输入为扭矩，输出为位置的积分过程进行频率分析。如果输入是正弦波，输出就变成这样：
我用来提取振幅比和相位的代码如下所示：
freq = 40;
freq_rad = freq * 2 * pi
phase_offset_rad = 30 * pi / 180
gain = 0
fs = 500;
L = 100;
t = (0:L-1)*(1/fs);
in = 2 * sin(freq * 2 * pi * t);
pos_in = [];
vel = 0;
pos = 0;
for i = 1:length(t)
vel = vel + in(i);
pos = pos + vel;
pos_in = [pos_in; pos];
end
out = pos_in;
%out = (2 + gain) * sin(freq * 2 * pi * t + phase_offset_rad);
fft_in = fft(in);
fft_out = fft(out);
[mag_in idx_in] = max(abs(fft_in));
[mag_out idx_out] = max(abs(fft_out));
phase = angle(fft_out(idx_out)) - angle(fft_in(idx_in))
phase_deg = phase / (pi / 180)
ratio = mag_out / mag_in
如果我 运行 它在完全直的正弦信号上那么它就可以完美地工作。但是，一旦我像上面那样添加失真，相位和幅度值就不对了。我想我需要以某种方式 "flatten" 信号。但我不确定如何从中提取正确的振幅。振幅是多少？我会在输出中说，从一个 "plateau" 到下一个 "plateau" 的测量值约为 45，因为这就是物体移动的距离。那将是〜22.5的比率。然而计算的结果是196。
可能是我想错了？我想最终使用实验数据导出从扭矩输入到位置输出的传递函数。也许有人可以展示如何做到这一点？
我一直在想，我能做的就是记录振幅比和相位，然后绘制波德图，从中轻松提取传递函数。到目前为止，我还无法从具有不同输入频率的 运行ning 测试中得到波德图。
因为 FFT 假设您对完全周期性信号（正好是信号的一个周期）执行频率分析，所以您的 fft(out) 将包含非常大的功率干扰（参见 Periodicity and Shift theorem）。
我相信，您的案例可以通过执行一些系统修改来避免 FFT 分析伪影。您可以估算系统+滤波器的传递函数，而不是估算系统的传递函数。 IE。您必须通过高通滤波器传递系统的输出信号：
out = filter([1 -1], 1, out);
然后，您就可以进行分析了。
滤波器的频率响应可以通过 freqz 函数来估计（或者在我的例子中只是 H = fft([1 -1], length(out));）。然后，您可以通过应用 fft_out = fft_out ./ H(:); 的逆响应来消除频域中的滤波器影响。另外，不要忘记在最大估计之前取消第 0 个频率 fft_out(1) = 0;。
顺便说一下，不同频率的相位差估计看起来很奇怪（在你的代码中phase = angle(fft_out(idx_out)) - angle(fft_in(idx_in))）。看起来你必须使用 idx_in（或 idx_out，取决于哪种估计更可靠）来计算回合角度。
注意：此答案并非完整指南，可能需要一些现实生活中的增强功能。
P.S。尝试在实际应用中应用开窗进行频率响应估计（例如 Hamming window）。
P.P.S 尝试在 https://dsp.stackexchange.com/
上提问
更新：
在某些情况下，您可以对输入信号执行相同的线性输入信号变换，而不是忽略滤波器的影响：in = filter([1 -1], 1, in);
我需要对一个输入为扭矩，输出为位置的积分过程进行频率分析。如果输入是正弦波，输出就变成这样：
我用来提取振幅比和相位的代码如下所示：
freq = 40;
freq_rad = freq * 2 * pi
phase_offset_rad = 30 * pi / 180
gain = 0
fs = 500;
L = 100;
t = (0:L-1)*(1/fs);
in = 2 * sin(freq * 2 * pi * t);
pos_in = [];
vel = 0;
pos = 0;
for i = 1:length(t)
vel = vel + in(i);
pos = pos + vel;
pos_in = [pos_in; pos];
end
out = pos_in;
%out = (2 + gain) * sin(freq * 2 * pi * t + phase_offset_rad);
fft_in = fft(in);
fft_out = fft(out);
[mag_in idx_in] = max(abs(fft_in));
[mag_out idx_out] = max(abs(fft_out));
phase = angle(fft_out(idx_out)) - angle(fft_in(idx_in))
phase_deg = phase / (pi / 180)
ratio = mag_out / mag_in
如果我 运行 它在完全直的正弦信号上那么它就可以完美地工作。但是，一旦我像上面那样添加失真，相位和幅度值就不对了。我想我需要以某种方式 "flatten" 信号。但我不确定如何从中提取正确的振幅。振幅是多少？我会在输出中说，从一个 "plateau" 到下一个 "plateau" 的测量值约为 45，因为这就是物体移动的距离。那将是〜22.5的比率。然而计算的结果是196。
可能是我想错了？我想最终使用实验数据导出从扭矩输入到位置输出的传递函数。也许有人可以展示如何做到这一点？
我一直在想，我能做的就是记录振幅比和相位，然后绘制波德图，从中轻松提取传递函数。到目前为止，我还无法从具有不同输入频率的 运行ning 测试中得到波德图。
因为 FFT 假设您对完全周期性信号（正好是信号的一个周期）执行频率分析，所以您的 fft(out) 将包含非常大的功率干扰（参见 Periodicity and Shift theorem）。
我相信，您的案例可以通过执行一些系统修改来避免 FFT 分析伪影。您可以估算系统+滤波器的传递函数，而不是估算系统的传递函数。 IE。您必须通过高通滤波器传递系统的输出信号：
out = filter([1 -1], 1, out);
然后，您就可以进行分析了。
滤波器的频率响应可以通过 freqz 函数来估计（或者在我的例子中只是 H = fft([1 -1], length(out));）。然后，您可以通过应用 fft_out = fft_out ./ H(:); 的逆响应来消除频域中的滤波器影响。另外，不要忘记在最大估计之前取消第 0 个频率 fft_out(1) = 0;。
顺便说一下，不同频率的相位差估计看起来很奇怪（在你的代码中phase = angle(fft_out(idx_out)) - angle(fft_in(idx_in))）。看起来你必须使用 idx_in（或 idx_out，取决于哪种估计更可靠）来计算回合角度。
注意：此答案并非完整指南，可能需要一些现实生活中的增强功能。
P.S。尝试在实际应用中应用开窗进行频率响应估计（例如 Hamming window）。
P.P.S 尝试在 https://dsp.stackexchange.com/
上提问更新： 在某些情况下，您可以对输入信号执行相同的线性输入信号变换，而不是忽略滤波器的影响：in = filter([1 -1], 1, in);