时间序列分析-提取并跟踪周期分量

Question

我的数据集的频谱显示了时间序列中的3个周期分量。我想减去周期性成分，以保持没有这些周期性的数据。它指出了周期性事件，其周期为（1/144=每天）、（1/72=1/2天）和（1/6=每小时）。我的想法是找出我的数据集的傅立叶分量（mag和相位），提取这3个特定频率的傅立叶分量，并创建一个新信号，其为：数据-PeriodicSignal_1h-PeriodicSignl_1/2天-Periodicsignl_1天我尝试使用fft，但我不知道如何提取这些特定频率的信号。

我的数据集很复杂，但我正在编写一个示例来理解这个过程。以下是示例：

samplingFrequency = 1000;
timeInterval = 1/samplingFrequency;
signalIndex = seq(0, 1, by=timeInterval);
N = 1000
a1 = 2;  
a2 = 3;  
f1 = 10; 
f2 = 20; 
signal1 = a1 * sin(2 * pi * f1 * signalIndex);
signal2 = a2 * sin(2 * pi * f2 * signalIndex);
inputSignal = signal1 + signal2;

Y<-fft（inputSignal）
mag&lt-sqrt（Re（Y）^2+Im（Y）^ 2）*2/长度（输入信号）
相位&lt-atan（Im（Y）/Re（Y））Yr<-Re（Y）Yi<-Im（Y）

我试图提取频率为f1的信号的mag和相位。我想生成一个新的信号是：

脉冲信号-Signal_f1

Answer 1

我相信以下内容符合您的要求。。。我更改了一些变量名。底部是您询问的频率选择。

设置时间和频率参数

samplingFrequency = 1000;
f_Hz = samplingFrequency
N = 1000
df_Hz = f_Hz / N
T = 1 / df_Hz
dt=T/N
t = dt*(seq(1,N)-1)

生成伪信号，无噪声

a1 = 2;  
a2 = 3;  
f1 = 10; 
f2 = 20; 
signal1 = a1 * sin(2 * pi * f1 * t);
signal2 = a2 * sin(2 * pi * f2 * t);
inputSignal = signal1 + signal2;

绘制伪信号

plot(t, signal1,type='l',col='green',ylim=c(-6,6))
lines(t, signal2,col='red')
lines(t, inputSignal,col='black')

获取fft，并绘制正频率部分

Y <- fft(inputSignal)
m <- floor(N/2)-1
posFreqIndices <- 2:(m+1)
negFreqIndices <- N:(m+3)
mag <-sqrt(Re(Y)^2+Im(Y)^2)*2/length(inputSignal)
phase <-atan(Im(Y)/Re(Y))
Yr <- Re(Y) 
Yi <- Im(Y)
freq <- seq(df_Hz,f_Hz/2-df_Hz,df_Hz)
plot(freq,mag[posFreqIndices],type='l',xlab='Freq (Hz)', ylab='Magnitude',xlim=c(0,30))
# plot(freq,10*log10(mag[posFreqIndices]),type='l',xlab='Freq (Hz)', ylab='Magnitude (db)',xlim=c(0,30))
# plot(freq,phase[posFreqIndices]*180/pi,type='l',xlab='Freq (Hz)', ylab='Phase (deg)',xlim=c(0,30))

根据幅度识别滤波信号的频率

ampSelectIndices <- which(mag>1.9 & mag < 2.1)

生成所选频率的滤波fft

YAmpSelect <- Y*0
YAmpSelect[ampSelectIndices] = Y[ampSelectIndices]

计算逆fft

yAmpSelect = Re(fft(YAmpSelect, inverse = TRUE))/length(YAmpSelect)

绘制滤波信号

plot(t,yAmpSelect,t='l',xlab='t (sec)',ylab='Filtered for mag ~ 2')

绘制原始信号减去滤波信号

plot(t,inputSignal-yAmpSelect,type='l')

fft是在频率折叠的情况下计算的。以下检查展开过程，此检查适用于实值信号（非复值时间信号）。该过程对于复值时间信号是正确的。

checkFreqWrapping = all.equal(mag[posFreqIndices], mag[negFreqIndices])
stopifnot(checkFreqWrapping)

通过频率选择fft值

freqSelectIndices_a <- which(9.95 < freq & freq < 10.05)
freqSelectIndices = union(posFreqIndices[freqSelectIndices_a],negFreqIndices[freqSelectIndices_a])

为选定频率创建fft

YFreqSelect <- Y*0
YFreqSelect[freqSelectIndices ] = Y[freqSelectIndices ]

计算时间信号，绘制它。

yFreqSelect = Re(fft(YFreqSelect, inverse = TRUE))/length(YFreqSelect)
plot(t,yFreqSelect,t='l',xlab='t (sec)',ylab='Filtered for mag ~ 2')
plot(t,inputSignal-yFreqSelect,type='l')

好吧，我想这解释了如何根据频率选择fft值。。。祝你好运