什么是 numpy.fft.rfft 和 numpy.fft.irfft 及其在 MATLAB 中的等效代码
What is numpy.fft.rfft and numpy.fft.irfft and its equivalent code in MATLAB
我正在将 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));
我误会了什么?
numpy 中的实数 FFT 使用实值函数的傅立叶变换 "skew-symmetric",即频率 k
处的值是该值的复共轭k=1..N-1
的频率 N-k
(正确的术语是 Hermitian)。因此 rfft
returns 只有结果中对应于非正频率的部分。
对于大小为 N
的输入,rfft
函数 returns FFT 输出的一部分对应于等于或低于 N/2
的频率。因此,如果 N
是偶数(从 0
到 N/2
的所有频率),rfft
的输出大小为 N/2+1
,如果 N
是奇数(从 0 到 (N-1)/2
的所有频率)。观察函数 floor(n/2+1)
returns 偶数和奇数输入大小的正确输出大小。
所以要在 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 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
我正在将 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));
我误会了什么?
numpy 中的实数 FFT 使用实值函数的傅立叶变换 "skew-symmetric",即频率 k
处的值是该值的复共轭k=1..N-1
的频率 N-k
(正确的术语是 Hermitian)。因此 rfft
returns 只有结果中对应于非正频率的部分。
对于大小为 N
的输入,rfft
函数 returns FFT 输出的一部分对应于等于或低于 N/2
的频率。因此,如果 N
是偶数(从 0
到 N/2
的所有频率),rfft
的输出大小为 N/2+1
,如果 N
是奇数(从 0 到 (N-1)/2
的所有频率)。观察函数 floor(n/2+1)
returns 偶数和奇数输入大小的正确输出大小。
所以要在 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 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