信息处理 - 生成具有任意功率谱的时间序列 - 吾爱随笔录

生成具有任意功率谱的时间序列

信息处理 fft 功率谱密度

2022-02-05 10:38:56

我需要使用给定的 PSD 生成随机数。为此，我找到了这个食谱。如果它在付费墙后面，现在它可以工作了：

给定您的功率谱与使用时间步长的时间序列生成的点，您定义 $P_k$ $0\leq k<N$ $h$

A_{k} = \sqrt{\frac{P_{k}}{2 N h}} \exp (i a_{k})

$A_k=\sqrt{\frac{P_k}{2Nh}}\exp{(ia_k)}$

其中中均匀随机选择。的傅里叶变换。这里的想法是，由于 PSD 丢弃了相位信息，我们只需为变换生成随机相位，以获得生成 PSD 的时间序列的可能实现。 $a_k$ $[0,2\pi)$ $A_k$

但是，为了确保所有得到的时间点都是真实的，作者给出了的条件，或者等价。 $A_{N-k}=A^*_k$ $a_{N-k}=-a_k$

然而，当我实现这个时，我没有得到真正的信号，所以我似乎做错了什么。我马上看到了一个问题：对于，FFT 假定一个周期信号，因此满足条件的唯一方法是。然后您可以按如下方式生成其余阶段（在 python 中）： $k=0$ $a_0=0$

def gen_phase(N):
    phase = np.zeros(N)
    randoms = np.random.random(size=N/2)*2*np.pi
    phase[1:N/2+1] = randoms
    phase[N/2+1:] = np.flipud(-randoms)
    return phase

另一个问题是不可能满足给定的条件是是偶数，据我所知。这实际上很痛苦，因为 fft 对于奇数来说非常慢（尤其是 ... 形式的） $N$ $2^n+1$

有人可以建议我可能做错了什么吗？我得到的信号看起来或多或少是幅度的正确离子项，但它具有大致相等的实部和虚部，所以有些地方是错误的。其余的代码太简单了，我认为是罪魁祸首：

def main():
    psd = read_psd()
    N = len(psd)
    phase = gen_phase(N)
    h = 0.24e-6
    A = np.sqrt(psd/(2*N*h)) * np.exp(1j*phase)
    timeseries = np.fft.fft(A)
    pl.plot(np.real(timeseries))
    pl.show()

if __name__=='__main__':
    main()

编辑：我误解了。因为取决于，所以我首先没有选择的自由，因为我无法控制幅度。有任何想法吗？ $A_k$ $P_k$ $A_{N-k}=A^*_k$

3个回答

希望这个 octave 实现可以帮助您：

N = 1e3;
P = ones(N,1);
Ts = 1/1000;
A = sqrt(P/(2*N*Ts)) .* exp( 1j * rand(N,1) );

% DC sample must be real.
n0 = floor( N/2 ) + 1;
A(n0) = P(n0);

% For even-length cases, the first sample has no
% equivalent positive frequency, so we force it
% to be real too.
A(1) = abs( A(1) );

% Force the spectrum to be conjugate symetric.
if( mod(N,2) == 0 )
  A(2:n0-1) = flipud( conj(A(n0+1:end)) );
else
  A(1:n0-1) = flipud( conj(A(n0+1:end)) );
end

% Time domain data.
x = ifft( ifftshift( A ) );

你正在颠倒你的阶段。如果将 N/2 的相位设置为零，您应该能够使用 N 的偶数值。

您遇到的问题是您的 PSD 跨越整个 N。除非否则您不会生成具有复杂共轭对称性的 DFT，而这是产生实值信号所需的。 $P_k = P_{N-k}$

此外，您应该使用“ifft”而不是“fft”来进行逆 FFT。

希望这可以帮助，

赛德

编辑：我在发布之前没有看到您的编辑。看来你想通了。

=======================

跟进：

对于偶数 N：

定义 gen_phase(N):
    相位 = np.zeros(N)
    随机数 = np.random.random(size=N/2-1)*2*np.pi
    相位[1:N/2] = 随机数
    相位[N/2+1:] = np.flipud(-randoms)
    返回阶段

感谢与其他回答者的讨论，我能够弄清楚。我在这里发布我自己的答案，因为其他人都没有完整的图片，尽管他们确实有帮助。

我缺少的技巧是，因为我的输入 PSD 仅针对正频率定义，所以我必须首先将 PSD 对称地复制到 0 左右。然后，由于 numpy 的 fft 实现内部频率分量的特定顺序，需要 ifftshift让一切都顺理成章。

这是完整的代码，假设您将 PSD 保存在为 [0, fmax] 中的频率定义的 .csv 文件中。

import numpy as np
from scipy.fftpack import fft, fftshift, ifftshift
from scipy.signal import welch
import pandas as pd
import matplotlib.pyplot as pl

def read_psd():
    psd = pd.read_csv('psd-example.csv',header=None, names=['f','P','I'])
    return psd['f'],psd['P']

def gen_phase(N):
    if N%2==0:
        raise ValueError('THis methofd requires a symmetric PSD, which can only be achieved with an odd number of data points')
    else:
        phase = np.zeros(N)
        randoms = np.random.random(size=N/2)*2*np.pi-np.pi
        phase[:N/2] = randoms
        phase[N/2+1:] = np.flipud(-randoms)
    return phase


def symmetrize_psd(psd):
    N = len(psd)
    symmetric = np.zeros(2*N-1)
    symmetric[0:N] = np.flipud(psd)
    symmetric[N:] = psd[1:]
    return symmetric

def scale_psd(psd, dt):
    N = len(psd)
    phase = gen_phase(N)
    A = np.sqrt(psd/(2*N*dt)) * np.exp(1j*phase)
    return ifftshift(A)

def get_timeseries(A):
    timeseries = np.fft.fft(A)
    return np.real(timeseries)

def main():
    averages = 200
    fi, psd = read_psd()
    psd = psd[:-1]
    fi=fi[:-1]

    sympsd = symmetrize_psd(psd)
    N = len(sympsd)
    timeseries = np.empty(N*averages)
    for i in range(averages):
        A = scale_psd(sympsd, 0.24e-6)
        timeseries[i*N:(i+1)*N] = get_timeseries(A)
    f, pxxr = welch(np.real(timeseries), fs=4166666.67, nperseg=2**np.floor(np.log2(N)+1))
    pl.loglog(f, pxxr)
    pl.loglog(fi, psd)
    pl.ylim((1e-5,100))
    pl.show()

if __name__=='__main__':
    main()

其它你可能感兴趣的问题

上一篇是否有任何用于自适应噪声消除的开源库？下一篇如何从 PCM 音频数据中获取音量级别？