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

j2datikz  于 2023-10-19  发布在  Matlab
关注(0)|答案(1)|浏览(177)

我正在将Python代码转换为MATLAB,其中一个代码使用了numpy rfft。在numpy的文档中,它说的是真实的输入。
计算真实的输入的一维离散傅立叶变换。
所以我在MATLAB中使用了abs,但结果不同。
Python代码

ffta = np.fft.rfft(a)

MATLAB代码

ffta = abs(fft(a));

我误解了什么?

hk8txs48

hk8txs481#

numpy中的真实的FFT使用了这样一个事实,即真实的值函数的傅立叶变换是“斜对称”的,也就是说,频率k处的值是频率N-k处的值对于k=1..N-1的复共轭(正确的项是 Hermitian)。因此rfft只返回对应于非正频率的结果部分。
对于大小为N的输入,rfft函数返回FFT输出中对应于N/2或以下频率的部分。因此,如果N是偶数(所有频率从0N/2),则rfft的输出大小为N/2+1,或者如果N是奇数(所有频率从0到(N-1)/2),则(N+1)/2。注意,函数floor(n/2+1)对于偶数和奇数输入大小都返回正确的输出大小。
在matlab中实现rfft

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 out = irfft(x, even)
    if nargin == 1
        even = true;
    end
    % `n`: the output length
    % `s`: the variable that will hold the index of the highest
    %      frequency below N/2, s = floor((n+1)/2)
    N = numel(x);  % function assumes 1D input
    if even
        n = 2 * (N - 1);
        s = N - 1;
    else
        n = 2 * (N - 1) + 1;
        s = N;
    end
    outf = zeros(1, n);
    outf(1:N) = x;
    outf(N+1:n) = conj(x(s:-1:2));
    out = ifft(outf);
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

相关问题