在 R 中重建 FFT 处理的噪声信号

Question

您好，我想隔离所有由 R 中的快速傅里叶变换 (FFT) 产生的噪声信号的正弦和余弦。这是为了说明 FFT 在小且定期采样的时间序列上的噪声信号行为。我从 Matlab 的解释中导出了一个脚本； https://nl.mathworks.com/help/matlab/ref/fft.html

通过简单的正弦相加，脚本表现良好：

# Create a signal with given parameters ----

L <- 1500  # Length of data

Fs <- 1000 # Sampling frequency

Ts <- 1/Fs # Sampling rate

t <- (0:(L-1))*Ts # Time value

S1 <- 0.7*sin(2*pi*50*t)
S2 <- sin(2*pi*120*t)

S <- S1 + S2

X <- S

# Uncomment to add noise ----
# set.seed(42)
# X <- S + 0.5*rnorm(length(t))

# Perform FFT on X ----

Y <- fft(X)

r1 <- Re(Y)/L
i1 <- Im(Y)/L

# Rearrange fft output to get the frequency, 
# the real and the imaginary components well identified ----

r1 <- r1[1:((L/2)+1)]
r1[2:(length(r1)-1)] <- 2 * r1[2:(length(r1)-1)]

i1 <- i1[1:((L/2)+1)]
i1[2:(length(i1)-1)] <- 2 * i1[2:(length(i1)-1)]

f <- Fs*(0:(L/2))/L

time <- t
freq <- f
real <- r1
imag <- i1

# Reconstitute each and every sine and cosine ----

lt <- length(time)
lf <- length(freq)

mtime <- matrix(rep(time, lf), nrow = lt)
mfreq <- matrix(rep(freq, lt), nrow = lt, byrow = T)

mcos <- cos(2 * pi * mtime * mfreq)
msin <- sin(2 * pi * mtime * mfreq)

acos <- matrix(rep(real, each = lt), nrow = lt)
asin <- matrix(rep(imag, each = lt), nrow = lt)

rcos <- mcos * - acos # Negative for whatever reason
rsin <- msin * - asin # Negative for whatever reason

# Add real and imaginary parts (cosines and sines) ----

comb <- rcos + rsin 

# Reconstitute the entire signal ----

synth <- rowSums(comb)

# Plot ----

par(mfrow = c(1,4))

ylim <- c(0,0.2)
xlim <- NULL

plot(X, time, type = "l", pch = 19, xlab = "Signal", 
     xlim = xlim, ylim = ylim)

# 181 index of f = 120
plot(comb[,181] ,time, type = "l", xlab = "Isolated frequencies", 
     xlim = xlim, ylim = ylim, lty = 5)

# 76 index of f = 50
lines(comb[,76] ,time, type = "l", lwd = 2)

plot(synth ,time, type = "l", xlab = "Reconstituted signal", 
     xlim = xlim, ylim = ylim)

difference <- synth - X

hist(difference, breaks = 100, col = "black")

这段代码给出了下图。左边的图是我应用 FFT 的信号，左中的一个是组成信号的两个正弦，由 FFT 提取，右中的图是所有正弦曲线的相加.信号和噪声之间的差异通过右侧的直方图来表征。它非常小，所以我假设这是浮点算术错误的结果，可以忽略不计。

我的问题是当我处理高噪声信号时； FFT 重建与初始信号明显不同，如下图所示（与之前相同的解释，相同的代码，我只是取消了添加噪声的代码位）。

重建的信号明显不同于初始信号（尽管具有明显相同的方差，并且添加了相同的正弦）。这个问题可以避免吗？

Answer 1

好的，它可以使用直接傅里叶变换（冰川傅里叶变换？）

# Create a signal with given parameters ----

L <- 1500  # Length of data

Fs <- 1000 # Sampling frequency

Ts <- 1/Fs # Sampling rate

dt <- (0:(L-1))*Ts # Time value

S1 <- 0.7*sin(2*pi*50*dt)
S2 <- sin(2*pi*120*dt)

S <- S1 + S2

xy <- S

# Uncomment to add noise ----
set.seed(42)
xy <- S + 0.5*rnorm(length(t)) + 20

# Glacial/Direct FOurier Transform

fseq <- Fs*(0:(L/2))/L

ahr <- NULL
ahi <- NULL

for(f in fseq){

  hr <- 2*sum(xy * cos(2 * pi * f * dt))/L
  hi <- 2*sum(xy * sin(2 * pi * f * dt))/L

  ahr <- c(ahr, hr)
  ahi <- c(ahi, hi)

}

ahr[1] <- ahr[1]/2
ahi[1] <- ahi[1]/2

ahr[length(ahr)] <- ahr[length(ahr)]/2
ahi[length(ahi)] <- ahi[length(ahi)]/2

time <- dt
freq <- fseq
real <- ahr
imag <- ahi

# Reconstitute each and every sine and cosine ----

lt <- length(time)
lf <- length(freq)

mtime <- matrix(rep(time, lf), nrow = lt)
mfreq <- matrix(rep(freq, lt), nrow = lt, byrow = T)

mcos <- cos(2 * pi * mtime * mfreq)
msin <- sin(2 * pi * mtime * mfreq)

acos <- matrix(rep(real, each = lt), nrow = lt)
asin <- matrix(rep(imag, each = lt), nrow = lt)

rcos <- mcos * acos
rsin <- msin * asin

# Add real and imaginary parts (cosines and sines) ----

comb <- rcos + rsin

# Reconstitute the entire signal ----

synth <- rowSums(comb)

# Plot ----

par(mfrow = c(1,4))

ylim <- c(0,0.2)

plot(xy, time, type = "l", pch = 19, xlab = "Signal",
     ylim = ylim)

# 181 index of f = 120
plot(comb[,181] ,time, type = "l", xlab = "Isolated frequencies",
     ylim = ylim, lty = 5)

# 76 index of f = 50
lines(comb[,76] ,time, type = "l", lwd = 2)

plot(synth ,time, type = "l", xlab = "Reconstituted signal",
     ylim = ylim)

difference <- synth - xy

hist(difference, breaks = 100, col = "black")