【发布时间】:2021-07-30 16:55:14
【问题描述】:
我正在尝试使用 iterative Radix2 算法 实现我自己的 FFT(我什至考虑过使用递归算法,但似乎 Simulink 不允许我这样做)。 我尝试了它的所有可能变体,比如 DIT 和 DIF,但我总是遇到同样的问题:输出对我来说没有多大意义。 它们具有相似的输出,包括 DIF 和 DIT,这两种都是用不同的代码变体实现的。
实际应用需要此算法,信号将被缓冲,然后将这些数据包窗口化并重叠。然后,窗口函数将被平均并发送到 Simulink 中的下一个模块,以获得 Welch 谱估计。
这是我实现的(其中一个)DIT 版本的代码:
function [FFT,Wk,k] = radix2(u,k_old,lap) %lap = 512 = 1024/2 with 1024 Window size
N = length(u);
P = log2(N);
FFT = zeros(N,1)+0*1j;
omega = 0j;
% Hann Window
n1 = (k_old)*lap:(N+k_old*lap-1);
w1 = sin(pi.*n1/(N-1)).^2;
z1 = u.*w1';
z1=z1';
% Average
Wk = 0;
for i=1:N
Wk = Wk + w1(i)^2;
end
%% Algorithm
A = z1(bitrevorder(1:N));
for s=1:log2(N)
m=2^s;
Wm = exp(2*pi*1j/m)
W = 1 +0j;
for j= 1:m/2
for k=j:m:(N)
odd = W*A(k+m/2)
even = A(k);
FFT(k)=odd + even;
FFT(k+m/2) =even-odd;
end
W = W*Wm;
end
end
这是输出
而不是有两个/三个不同的尖峰(输入是带有一些噪声的正弦信号,即使我关闭噪声或窗口函数,结果也没有太大变化,所以这些是正确的),我得到这个蝴蝶形的信号。为什么会这样?
我一开始以为是频率归一化问题,但是没有解决。
我真的希望任何人都可以提供帮助,这个问题不是太愚蠢,我现在完全迷路了。
重新编辑:
原因是我在动态参数估计的不同技术之间进行并行处理:锁定放大器、Welch with DFT、FFT、Goertzel 和其他方法。然后我会比较方差、时间、FLOPS 之类的东西......所以自我实现它正是我工作的重点。我知道 welch 已经存在于 Matlab 和 Simulink 中,但是对于我的特殊研究目的,它没有给我带来任何好处。此外,我的大多数版本都运行良好,唯一的一个问题是 迭代 FFT 实现,它看起来不像傅立叶变换,我真的不明白为什么。
输出应在 + 和 -50Hz 处显示两个尖峰,一个在 0Hz 处,因为信号的比例项和其他一些较小的周围,与 +-50Hz 和 0 相比微不足道。
【问题讨论】:
-
为什么不在 MATLAB 中使用
fft函数?为什么要重新实现它?人们制作自己的实现来了解它是如何工作的,但不是为了在实际应用程序中使用。使用您的软件随附的非常高效、经过良好测试的实现。 -
正如我在评论中所说的,我想研究可能的实现及其结构,计算失败次数等等,所以我确实需要整个实现,如果它更糟也没关系.它必须适合比较论文,并很好地解释为什么它比另一个更好或更差。当问题是纯粹的结构性问题时,为什么知道我为什么需要它如此重要?你看,即使 radix2 实现的 Flops 通常被计算为 N*log2(N),结构本身也可以稍微优化或不优化算法,我想看看这个。还是谢谢你。
-
如果我正在审阅你的论文,我也会说同样的话。你在比较 A 和 B,但是 B 使用了一个糟糕的自定义 FFT 实现,那么比较是不公平的。在多线程、SIMD 指令、流水线等时代,FLOPS 变得无关紧要。使用优化的实施并比较时间、能源使用或任何相关的内容。
-
另外,我不是在问你需要它做什么,这是一个反问。我告诉你不要这样做。这是一个来自几十年经验的建议。不要重新发明轮子,也不要重新实现 FFT。
标签: matlab fft simulink spectrum radix-sort