fft() 的 Matlab 和 R 函数：N 点 FFT 和维度问题

信息处理 fft 频谱自习

2022-02-21 01:50:23

在 MATLAB 中，Y = fft(X,n,dim)返回沿维度的傅里叶变换dim。例如，如果X是一个矩阵，则fft(X,n,2)返回每一行的 n 点傅里叶变换。但是，在 R 中是不同的。你知道解决这个问题的任何功能吗？所以考虑n和dim，得到一个n长度复向量？请参阅示例中的那个n>ncol(X).

MATLAB代码

Fs = 1000; % Sampling frequency
T = 1/Fs; % Sampling period
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector

x1 = cos(2*pi*50*t); % First row wave
x2 = cos(2*pi*150*t); % Second row wave
x3 = cos(2*pi*300*t);  % Third row wave

n = 2^nextpow2(L);
dim = 2;
X = [x1; x2; x3];

Y_right = fft(X,n,dim); %Result: 3*1024 complex double

 Y_wrong = fft(X); %Result: 3*1000 complex double (I am looking for the above result)

R代码

Fs <- 1000
T <- 1/Fs
L <- 1000
t <- (0:(L-1))*T

x1 <- cos(2*pi*50*t)
x2 <- cos(2*pi*150*t)
x3 <- cos(2*pi*300*t)

n<-2^ceiling(log2(abs(L)))
dim <- 2
X <- t(matrix(c(x1, x2, x3), ncol=3))

Y <- fft(X) 
# Result: 3*1000 cplx and it is not the same as `Y_wrong` (operates along each row) 
#To get Y_wrong:

Y <- matrix(0, nrow=3, ncol=1000)

for (i in 1:ncol(X)) {
Y[,i] <- fft(X[,i], inverse = T)
}

我的目标是得到Y_right，和Y认为通过使用apply我们得到的一样Y_wrong。先感谢您。

2个回答

R 允许您通过使用apply系列函数来循环遍历行或列。请注意，apply不一定比循环快，但更符合 R 的函数式编程性质。还有其他apply函数，例如lapply、vapply、sapply等，但出于您的目的，apply就足够了。另外，请注意，您可能必须使用t(result). 1 表示遍历行，2 表示遍历列。所以，如果 X 是矩阵，

apply(X, 1, fft) does the fft of the rows.

apply(X, 2, fft) does the fft of the columns.

请注意，使用 .row 或 .col 是我的惯例。我本可以使用 x、y 或 z，因为从 apply 的角度来看，这只是一个虚拟变量。

此外，您可以输入? applyR 以获得应用帮助。有时它们可能有点晦涩难懂，但至少它们会告诉您可用的论点。

我终于得到了解决方案。通过执行（matlab）查看函数代码，open ftt我们可以看到：

%   FFT(X,N) is the N-point FFT, padded with zeros if X has less
%   than N points and truncated if it has more.
%
%   FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the
%   dimension DIM.

因此，如果X少于n点，则向量应在末尾用零填充：

padding_X<- matrix(0, ncol=n-L, nrow=3)  
X<-cbind(X,padding_X)

然后我们得到一个 3*1024 的矩阵，用零填充。下面要做的是fft：

Y <- matrix(0, nrow=3, ncol=ncol(X))

for (i in 1:nrow(X)) {
  Y[i,] <- fft(X[i,], inverse = T)
}

而且，瞧：

其它你可能感兴趣的问题

上一篇如何选择低通巴特沃斯滤波器的阶数和截止频率？下一篇从亮度值恢复 RGB 值？