如何在没有相位跳跃的情况下连接正弦波

Question

我需要制作一个 python 脚本，它生成给定频率的正弦波并使用 pyaudio（阻塞模式）播放它们，我还需要能够在它运行时更改该频率，对其进行调制并使用 pyqtgraph 绘制它。现在我有一个生成数据块的线程，我对 'connect' 这些正弦的方法是获取 fft 然后计算角度 (numpy.angle)，将其存储在一个变量中并将其用作相位偏移到下一个块，但我没有得到预期的结果，也许我遗漏了一些东西或将它们混淆了。

import matplotlib.pyplot as plt
import numpy as np
import pyaudio

#-----------------------
CHUNK = 1024
RATE =  44100
CHANNELS = 2
FORMAT = pyaudio.paFloat32
#-----------------------

samples = int(CHUNK)
t = np.arange(samples) / RATE
con = 0


def generate_sine(a: float = 0.5, freq: float = 440.0):

    global con

    sine = a * np.sin(2.0 * np.pi * freq * t + con)

    # get the angle of the wave

    phase = np.angle(np.fft.fft(sine))

    # update ref var to generate subsequent sines
    # begining where the last ended

    con = phase[-1]

    return sine


def play_sine(data):

    pa = pyaudio.PyAudio()

    stream = pa.open(format=FORMAT,
                         channels=CHANNELS,
                         rate=RATE,
                         input=False,
                         output=True,
                         frames_per_buffer=CHUNK)

    stream.write(np.array(data).astype(np.float32).tostring())

    stream.close()

if __name__ == '__main__':

    f = 80

    chunks = generate_sine(freq=f)

    for i in range(0,4):

        chunks = np.concatenate((chunks, generate_sine(freq=f)))

    #for i in range(0,10):

    #play_sine(chunks)

    plt.plot(chunks)

    plt.show()

演示剧情

您可以在链接的图像中看到 x=1024、x=2048 等处存在不连续点。

Answer 1

以下解决方案并非 100% 令人满意，但如果您能够在关节处丢弃每个片段的一些信号，则似乎效果很好。

它使用Hilbert transform重建相位。不幸的是，在信号的边缘附近往往存在伪影（例如，请参见图中红色信号的最末端），我只能通过切断它们来消除这些伪影。

>>> import numpy as np
>>> from scipy.signal import hilbert
>>> 
# create two sinewaves and join them discarding some signal near the joint
# this is the black curve in the plot
>>> A = np.sin(np.linspace(0, 30, 10000))
>>> B = np.cos(np.linspace(0, 15, 10000))
>>> Cjump = np.concatenate([A[:-1000], B[1000:]])
>>> 
# join the same sinewaves using Hilbert
>>> AP = np.unwrap(np.angle(hilbert(A)))
>>> BP = np.unwrap(np.angle(hilbert(B)))
>>> CP = np.concatenate([AP[:-1000], BP[1000:] - BP[1000] + AP[-1000]])
>>> C = np.cos(CP)
# this is the red curve in the plot

为什么要进行希尔伯特变换？因为它可以在某种程度上处理调制信号：示例（蓝色原图，重建阶段的绿色）：

生成信号的代码：

# original (blue):
>>> A = np.sin(np.convolve(np.random.uniform(0.2, 1, 10).repeat(1000), np.ones(10)/1000, 'same').cumsum())
>>> AP = np.unwrap(np.angle(hilbert(A)))
# recon (green):
>>> AR = np.cos(AP)

Answer 2

您正在使用

生成信号

    a * sin(2πf * t + con)

其中 t 范围超过 [0 .. CHUNK/RATE)。

当您开始下一个区块时，t 将重置为零。要生成连续波形，您需要修改 con 以生成与前一个样本相同的结果相位值。

使用 FFT 是行不通的，因为您生成的信号不是样本的精确倍数 window，而且您实际上对样本末尾的相位感兴趣 window, 不是开始的阶段。

相反，您只需要 生成函数的 在 t=t_end 处的相位值，模 2π。

也就是说，您可以简单地使用：

con = 2.0 * np.pi * f * CHUNK/RATE + con

但该值会增长，如果您将许多块连接在一起，频率很高，则可能会导致最终的数值问题。由于正弦函数是周期性的，您只需将结束相位归一化到 0 到 2π 范围内：

con = math.fmod(2.0 * np.pi * f * CHUNK/RATE + con, 2.0 * np.pi)

如果将生成函数修改为：

    a * sin(2π * (f * t + con))

则con表示从上一个卡盘结转的一个完整周期的分数，可以避免除以2π的模数，这可能会稍微提高精度。

con = math.modf(f * CHUNK/RATE + con)[0]

尝试更清晰的解释：

注意：此技术之所以有效，是因为您完全知道前一个区块的生成方程式，并且正在生成下一个区块。如果这些条件中的任何一个发生变化，您将需要一种不同的技术来匹配块。

前一个块是使用以下序列的 sin() 生成的：

2πf₁*0/RATE+con₁, 2πf₁*1/RATE+con₁, ..., 2πf₁*1022/RATE+con₁, 2πf₁*1023/RATE+con₁

应该清楚的是，为了顺利过渡到下一个块，该块应该从 2πf₁*1024/RATE+con₁

的 sin() 开始

下一个块从 2πf₂*0/RATE+con₂ 的 sin() 开始。

因此，如果满足以下条件，我们将顺利过渡：

2πf₂*0/RATE + con₂ = 2πf₁*1024/RATE + con₁

或

     0      + con₂ = 2πf₁*1024/RATE + con₁

或

              con₂ = 2πf₁*1024/RATE + con₁

可以在您的 generate_sine 函数中写成：

con = 2.0 * np.pi * f * CHUNK/RATE + con

这是我上面回答中 "out of nowhere" 的等式。从那时起，由于 sin() 函数是 2π 周期性的，我只是执行模 2π 约简，以防止 sin() 的参数无限增长，导致数值不准确。

希望这能让事情变得更清楚。

如何在没有相位跳跃的情况下连接正弦波

How to concatenate sine waves without phase jumps

python

trigonometry

signal-processing

numpy

pyaudio