Как объединить синусоиды без фазовых скачков

Question

Мне нужно сделать скрипт на python, который генерирует синусоидальные волны заданной частоты и воспроизводит их, используя pyaudio (режим блокировки), мне также нужно иметь возможность изменять эту частоту во время работы, модулировать ее и отображать на графике с помощью pyqtgraph. На данный момент у меня есть поток, генерирующий порции данных, и мой подход к «соединению» этих синусов заключался в том, чтобы получить 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

Вы генерируете свой сигнал с

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

, где t находится в диапазоне [0 .. CHUNK/RATE).

Когда вы запускаете следующий блок, t сбрасывается в ноль.Чтобы создать непрерывный сигнал, вам нужно изменить con, чтобы сгенерировать то же значение результирующей фазы, что и в предыдущем примере.

Использование БПФ не будет работать, поскольку генерируемый вами сигнал не являетсяточное значение, кратное окну сэмпла, плюс вам действительно интересна фаза в конце окна сэмпла, а не фаза в начале.

Вместо этого вам просто нужна функция генерации значение фазы при 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]

---

Попытка aБолее четкое объяснение:

Примечание. Этот метод работает только потому, что вы точно знаете уравнение генерации для предыдущего фрагмента и генерируете следующий фрагмент.Если какое-либо из этих условий изменится, вам понадобится другая техника для соответствия кускам.

Предыдущий блок генерируется с использованием sin() следующей последовательности:

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

Itдолжно быть ясно, что для плавного перехода к следующему чанку этот чанк должен начинаться с sin() из 2πf₁*1024/RATE+con₁

Следующий чанк начинается с sin() из 2πf₂*0/RATE+con₂.

Поэтому у нас будет плавный переход, если:

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

- это уравнение, которое появилось "из ниоткуда" в моем ответе выше.С этой точки зрения, поскольку функция sin() является периодической 2π, я просто выполняю уменьшение по модулю 2π, чтобы аргумент sin() не увеличивался без границ, вызывая неточности чисел.

Надеюсь, что это делает вещияснее.

Answer 2

Следующее решение не является удовлетворительным на 100%, но, кажется, работает хорошо, если вы можете отбросить немного сигнала каждого фрагмента в суставах.

Используется преобразование Гильберта для восстановления фазы. К сожалению, вблизи краев сигнала, как правило, присутствуют артефакты (см., Например, самый конец красного сигнала на графике), от которого я мог избавиться только путем их обрезания.

>>> 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)