python - 什么是 numpy.fft.rfft 和 numpy.fft.irfft 及其在 MATLAB 中的等效代码

Question

我正在将 python 代码转换为 MATLAB，其中一个代码使用 numpy rfft。在 numpy 的文档中，它说的是真实输入。

Compute the one-dimensional discrete Fourier Transform for real input.

所以我在 MATLAB 中所做的是使用 abs 但结果不同。

Python代码

ffta = np.fft.rfft(a) 

MATLAB 代码

ffta = abs(fft(a));

我误解了什么？

Answer 1

numpy 中的真实 FFT 使用这样一个事实，即实值函数的傅立叶变换可以说是“偏斜对称”，即频率 k 处的值。是频率 N-k 处的值的复共轭为 k=1..N-1 (正确的术语是 Hermitian)。因此rfft只返回结果中对应于非正频率的部分。

对于大小为 N 的输入rfft函数返回对应于等于或低于 N/2 的频率的 FFT 输出部分.因此 rfft 的输出尺寸为 N/2+1如果 N是偶数(从 0 到 N/2 的所有频率)，或 (N+1)/2如果 N是奇数(从 0 到 (N-1)/2 的所有频率)。观察函数floor(n/2+1)为偶数和奇数输入大小返回正确的输出大小。

所以要重现rfft在matlab中

function rfft = rfft(a)
     ffta = fft(a);
     rfft = ffta(1:(floor(length(ffta)/2)+1));
end

例如

a = [1,1,1,1,-1,-1,-1,-1];
rffta = rfft(a)

会产生

rffta =

 Columns 1 through 3:

   0.00000 + 0.00000i   2.00000 - 4.82843i   0.00000 + 0.00000i   

 Columns 4 through 5:

   2.00000 - 0.82843i   0.00000 + 0.00000i

现在将其与 python 进行比较

>>> np.fft.rfft(a)
array([ 0.+0.j        ,  2.-4.82842712j,  0.-0.j        ,  
        2.-0.82842712j,  0.+0.j        ])

复制irfft

重现 irfft 的基本功能您需要从 rfft 恢复丢失的频率输出。如果所需的输出长度为偶数，则可以从输入长度计算输出长度为 2 (m - 1) .否则应该是 2 (m - 1) + 1 .

以下代码将起作用。

function irfft = irfft(x,even=true)
     n = 0; % the output length
     s = 0; % the variable that will hold the index of the highest
            % frequency below N/2, s = floor((n+1)/2)
     if (even)
        n = 2 * (length(x) - 1 );
        s = length(x) - 1;
     else
        n = 2 * (length(x) - 1 )+1;
        s = length(x);
     endif
     xn = zeros(1,n);
     xn(1:length(x)) = x;
     xn(length(x)+1:n) = conj(x(s:-1:2));
     irfft  = ifft(xn);
end

现在你应该有

>> irfft(rfft(a))
ans =

   1.00000   1.00000   1.00000   1.00000  -1.00000  -1.00000  -1.00000  -1.00000

并且

abs( irfft(rfft(a)) - a ) < 1e-15

对于奇数输出长度，你得到

>> irfft(rfft(a(1:7)),even=false)
ans =

   1.0000   1.0000   1.0000   1.0000  -1.0000  -1.0000  -1.0000