Как мне сгенерировать спектрограмму 1D сигнала в питоне?

Question

Я не уверен, как это сделать, и мне дали пример, например, спектрограмму, но это в 2D.

У меня есть код, который генерирует сочетание частот, и я могу выбрать их в FFT, как я могу увидеть их на спектрограмме? Я ценю, что частоты в моем примере не меняются со временем; значит ли это, что я увижу прямую линию через спектрограмму?

мой код и выходное изображение:

# create a wave with 1Mhz and 0.5Mhz frequencies
dt = 2e-9
t = np.arange(0, 10e-6, dt)
y = np.cos(2 * pi * 1e6 * t) + (np.cos(2 * pi * 2e6 *t) * np.cos(2 * pi * 2e6 * t))
y *= np.hanning(len(y))
yy = np.concatenate((y, ([0] * 10 * len(y))))

# FFT of this
Fs = 1 / dt  # sampling rate, Fs = 500MHz = 1/2ns
n = len(yy)  # length of the signal
k = np.arange(n)
T = n / Fs
frq = k / T  # two sides frequency range
frq = frq[range(n / 2)]  # one side frequency range
Y = fft(yy) / n  # fft computing and normalization
Y = Y[range(n / 2)] / max(Y[range(n / 2)])

# plotting the data
subplot(3, 1, 1)
plot(t * 1e3, y, 'r')
xlabel('Time (micro seconds)')
ylabel('Amplitude')
grid()

# plotting the spectrum
subplot(3, 1, 2)
plot(frq[0:600], abs(Y[0:600]), 'k')
xlabel('Freq (Hz)')
ylabel('|Y(freq)|')
grid()

# plotting the specgram
subplot(3, 1, 3)
Pxx, freqs, bins, im = specgram(y, NFFT=512, Fs=Fs, noverlap=10)
show()

Answer 1

То, что у вас есть, технически правильно, но вам просто нужно посмотреть на сигнал с интересной спектрограммой. Для этого вам нужно, чтобы частота менялась со временем. (И для того, чтобы это произошло, вам нужно много колебаний, поскольку требуется несколько колебаний, чтобы установить частоту, а затем вам нужно много из них, чтобы частота менялась со временем интересным образом.)

Ниже я изменил ваш код как можно меньше, чтобы получить частоту, которая делает что-то интересное (fscale просто нарастает частота со временем). Я публикую весь код, чтобы заставить его работать, но я изменяю только три из четырех верхних строк.

# create a wave with 1Mhz and 0.5Mhz frequencies
dt = 40e-9
t = np.arange(0, 1000e-6, dt)
fscale = t/max(t)
y = np.cos(2 * pi * 1e6 * t*fscale) + (np.cos(2 * pi * 2e6 *t*fscale) * np.cos(2 * pi * 2e6 * t*fscale))
y *= np.hanning(len(y))
yy = np.concatenate((y, ([0] * 10 * len(y))))

# FFT of this
Fs = 1 / dt  # sampling rate, Fs = 500MHz = 1/2ns
n = len(yy)  # length of the signal
k = np.arange(n)
T = n / Fs
frq = k / T  # two sides frequency range
frq = frq[range(n / 2)]  # one side frequency range
Y = fft(yy) / n  # fft computing and normalization
Y = Y[range(n / 2)] / max(Y[range(n / 2)])

# plotting the data
subplot(3, 1, 1)
plot(t * 1e3, y, 'r')
xlabel('Time (micro seconds)')
ylabel('Amplitude')
grid()

# plotting the spectrum
subplot(3, 1, 2)
plot(frq[0:600], abs(Y[0:600]), 'k')
xlabel('Freq (Hz)')
ylabel('|Y(freq)|')
grid()

# plotting the specgram
subplot(3, 1, 3)
Pxx, freqs, bins, im = specgram(y, NFFT=512, Fs=Fs, noverlap=10)
show()

Также обратите внимание, что здесь полезна только спектрограмма. Если вы видите хорошую форму волны или спектры, спектрограмма, вероятно, не будет интересна: 1) если форма волны ясна, у вас, вероятно, недостаточно данных и времени, в течение которого частота и четко определена, и изменяется достаточно, чтобы быть интересной; 2) если полные спектры ясны, у вас, вероятно, недостаточно частотная вариация для спектрограммы, так как спектр - это просто среднее значение того, что вы видите, изменяясь со временем на спектрограмме.

Если вы действительно хотите увидеть спектрограмму вашего исходного сигнала, вам просто нужно увеличить ось Y, чтобы увидеть ожидаемые пики (обратите внимание, что ось Y спектрограммы равна 2,5e8, должна быть больше, чем в вашем спектре):

Answer 2

Чтобы получить то, что вы ищете:

1) сэмплируйте сигнал 1d на высокой частоте (по крайней мере, в 5 раз превышающей частоту его наивысшей частотной составляющей)

2) использовать блоки выборок (степени 2, например, 1024,16384 и т. Д.) Для вычисления БПФ

3) для каждого спектра построить вертикальную линию пикселей, цвет которой представляет амплитуду каждой частоты.

4) повторите шаги 2 и 3 для каждого блока образцов.

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

Answer 3

Я только начинаю на Python 3.6 Спасибо за пример кода Spectrogram!

Однако, с Python 3.6 я немного изо всех сил пытался заставить этот пример кода спектрограммы работать (вызовы функций и деление чисел с плавающей запятой), я отредактировал код, так что теперь он работает на python 3.6 для моих приятелей-новичков на python.

наслаждаться

'''
Original Script for Python 2.7
https://stackru.com/questions/19052324/how-do-i-generate-a-spectrogram-of-a-1d-signal-in-python
Modified in August 2017 for Python 3.6
Python 2.7 two integers / Division generate Integer
Python 3.6 two integers / Division generate Float
Python 3.6 two integers // Division generate integer
'''


import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt


dt = 40e-9
t = np.arange(0, 1000e-6, dt)
fscale = t/max(t)
y = np.cos(2 * np.pi * 1e6 * t*fscale) + (np.cos(2 * np.pi * 2e6 *t*fscale) * np.cos(2 * np.pi * 2e6 * t*fscale))
y *= np.hanning(len(y))
yy = np.concatenate((y, ([0] * 10 * len(y))))

# FFT of this
Fs = 1 / dt  # sampling rate, Fs = 500MHz = 1/2ns
n = len(yy)  # length of the signal
k = np.arange(n)
T = n / Fs
frq = k / T  # two sides frequency range
frq = frq[range(n // 2)]  # one side frequency range
Y = fftpack.fft(yy) / n  # fft computing and normalization
Y = Y[range(n // 2)] / max(Y[range(n // 2)])

# plotting the data
plt.subplot(3, 1, 1)
plt.plot(t * 1e3, y, 'r')
plt.xlabel('Time (micro seconds)')
plt.ylabel('Amplitude')
plt.grid()

# plotting the spectrum
plt.subplot(3, 1, 2)
plt.plot(frq[0:600], abs(Y[0:600]), 'k')
plt.xlabel('Freq (Hz)')
plt.ylabel('|Y(freq)|')
plt.grid()

# plotting the specgram
plt.subplot(3, 1, 3)
Pxx, freqs, bins, im = plt.specgram(y, NFFT=512, Fs=Fs, noverlap=10)
plt.show()