我正在尝试实现来自DSP Guide的“合成方程”,方程8-2,这样我就可以在频域中对时间序列进行重采样。在我读方程的方法中,N是输出点的数量,给定k从0到N/2的循环,我最多只能重新采样到原始采样率的两倍。
我试着用R写了一个快速实现,但结果与我期望的一点也不接近。我的代码:
input <- c(1:9)
nin <- 9
nout <- 17
b <-fft(input)
reals <- Re(b) / (nout / 2)
imags <- Im(b) / (nout / 2)
reals[1] <- reals[1] / 2
reals[(nout/2)] <- reals[(nout/2)] / 2
output <- c(1:nout)
for (i in 1:nout)
{
realSum <- 0
imagSum <- 0
for (k in 1:(nout/2))
{
angle <- 2 * pi * (k-1) * (i-1) / nout
realSum <- realSum + (reals[k] * cos(angle))
imagSum <- imagSum - (imags[k] * sin(angle))
}
output[i] <- (realSum + imagSum)
}
对于我的输入(以1秒采样,重采样到0.5秒)
[1] 1 2 3 4 5 6 7 8 9
我得到了输出
[1] -0.7941176 1.5150954 0.7462716 1.5022387 1.6478971 1.8357487
2.4029773 2.1965426 3.1585254 2.6178195 3.7284660 3.3721128 3.8433588 4.6390705
[15] 3.4699088 6.3005605 2.8175240
而我的预期输出是
[1] 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9
我做错了什么?
发布于 2016-12-09 08:41:55
上面的代码看起来确实像您在文章中引用的逆向DFT。从您期望的输入/输出来看,您似乎需要一个廉价的线性插值,而不是DFT /逆DFT解决方案。
从上面所需的输入/输出关系语句可以看出,时间序列输入和所需时间序列输出之间具有所需的关系。
你编码的逆DFT将从频域到时域。
根据我从你的陈述中收集到的信息,从系列1,1->9开始,进行正向DFT (以获得频域),然后将结果进行反向DFT (以返回时域)。
希望这对你是一个正确的方向的提示。
https://stackoverflow.com/questions/41026273
复制相似问题