Здравствуйте, я хотел бы выделить все синусы и косинусы, которые создаются быстрым преобразованием Фурье (FFT) в R для шумовых сигналов. Это должно проиллюстрировать поведение БПФ с зашумленными сигналами на небольших и регулярно дискретизированных временных рядах. Я получил сценарий из объяснения 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")
Этот код дает следующий рисунок. График слева - это сигнал, к которому я применяю БПФ, в середине слева - два синуса, составляющих сигнал, извлеченный БПФ, а график в середине справа - это сложение всех синусоид. , Разница между сигналом и шумом характеризуется гистограммой справа. Он очень маленький, поэтому я предполагаю, что это результаты арифметических ошибок с плавающей запятой c ошибок и ничтожно малы.
Моя проблема - когда я работаю по сильно шумному сигналу; восстановление БПФ становится явно отличным от исходного сигнала, как можно видеть на следующем рисунке (то же самое объяснение, что и раньше, тот же код, я только раскомментировал бит кода, который добавляет шум).
Восстановленный сигнал явно отличается от исходного сигнала (несмотря на наличие, по-видимому, той же дисперсии и тех же синусов, которые были добавлены). Можно ли избежать этой проблемы?