Обнаружение пиков на графике БПФ

Мне было интересно, как можно обнаружить новые пики на графике БПФ в Python. скажем, у меня есть этот простой сюжет: введите описание изображения здесь И я хочу автоматически измерить «Сходство» или расположение пиков в зашумленном сигнале, я пытался использовать косинусное сходство, но мой реальный сигнал слишком шумно, и даже если я добавлю к сигналу новый пик, я продолжаю получать косинус 0,9, поскольку это только один пик. Это пример моего реального сигнала, и у меня также есть проблема, заключающаяся в том, что мой сигнал может смещаться в пределах измерений, поэтому я не могу получить стабильный массив частот, они могут находиться в окне +/- 100 Гц: enter Это код, который использовался для первого сюжета:

import numpy as np
from pylab import *
import scipy.fftpack

# Number of samplepoints
N = 600
# sample spacing
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
y1 = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)+ 0.7*np.sin(30.0 * 2.0*np.pi*x)+ 0.5*np.sin(10.0 * 2.0*np.pi*x)
y2 = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)+ 0.2*np.sin(60.0 * 2.0*np.pi*x)+ 0.4*np.sin(40.0 * 2.0*np.pi*x)
yf1 = scipy.fftpack.fft(y1)
yf2 = scipy.fftpack.fft(y2)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)

fig, ax = plt.subplots()
plot(xf, 2.0/N * np.abs(yf1[:N/2]))
plot(xf, 2.0/N * np.abs(yf2[:N/2]))
xlabel('Freq (Hz)',fontsize=16,weight='bold')
ylabel('|Y(freq)|',fontsize=16,weight='bold')
ax = gca()
fontsize = 14
for tick in ax.xaxis.get_major_ticks():
    tick.label1.set_fontsize(fontsize)
    tick.label1.set_fontweight('bold')
for tick in ax.yaxis.get_major_ticks():
    tick.label1.set_fontsize(fontsize)
    tick.label1.set_fontweight('bold')
grid(True)
show()
def cosine_similarity(v1,v2):
    "compute cosine similarity of v1 to v2: (v1 dot v2)/{||v1||*||v2||)"
    sumxx, sumxy, sumyy = 0, 0, 0
    for i in range(len(v1)):
        x = v1[i]; y = v2[i]
        sumxx += x*x
        sumyy += y*y
        sumxy += x*y
    return sumxy/math.sqrt(sumxx*sumyy)

print 'Cosine Similarity', cosine_similarity(2.0/N * np.abs(yf1[:N/2]),2.0/N * np.abs(yf2[:N/2]))

Я также хотел установить порог, но иногда пики в реальном сигнале могут быть меньше, чем предварительно определенный порог. Есть идеи ?


python python-2.7 plot similarity fft
person Mehdi    schedule 09.06.2016    source источник

Ответы (1)


arrow_upward
6
arrow_downward

Есть много способов найти пики и даже интерполировать их местоположение в подвыборке. Когда у вас есть пики, просто проверьте, не нашли ли вы новый.

Вы можете использовать пакет peakutils для поиска пиков. Там можно установить порог и минимальное расстояние между пиками.

import numpy as np
from pylab import *
import scipy.fftpack

# Number of samplepoints
N = 600
# sample spacing
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
y1 = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)+ 0.7*np.sin(30.0 * 2.0*np.pi*x)+ 0.5*np.sin(10.0 * 2.0*np.pi*x)
y2 = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)+ 0.2*np.sin(60.0 * 2.0*np.pi*x)+ 0.4*np.sin(40.0 * 2.0*np.pi*x)
yf1 = scipy.fftpack.fft(y1)
yf2 = scipy.fftpack.fft(y2)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)

v1 = 2.0/N * np.abs(yf1[:N/2])
v2 = 2.0/N * np.abs(yf2[:N/2])

# Find peaks
import peakutils
peaks_ind1 = peakutils.indexes(v1, thres=0.2, min_dist=5)
peaks_ind2 = peakutils.indexes(v2, thres=0.2, min_dist=5)

dist_th_for_new_peaks = 3
new_peaks = []
for p in peaks_ind2:
    found_new_peak = np.all(np.abs(p - peaks_ind1) > dist_th_for_new_peaks)
    if found_new_peak:
        new_peaks.append(p)
        print("New Peak!! - %d" % p)

fig, ax = plt.subplots()
plot(xf, v1, color='blue')
plot(xf, v2, color='green')
for p in peaks_ind1:
    ax.scatter(xf[p], v1[p], s=40, marker='s', color='blue', label='v1')
for p in peaks_ind2:
    ax.scatter(xf[p], v2[p], s=40, marker='s', color='green', label='v2')    
for p in new_peaks:
    ax.scatter(xf[p], v2[p], s=40, marker='s', color='red', label='new peaks')        

xlabel('Freq (Hz)',fontsize=16,weight='bold')
ylabel('|Y(freq)|',fontsize=16,weight='bold')

ax = gca()
fontsize = 14
for tick in ax.xaxis.get_major_ticks():
    tick.label1.set_fontsize(fontsize)
    tick.label1.set_fontweight('bold')
for tick in ax.yaxis.get_major_ticks():
    tick.label1.set_fontsize(fontsize)
    tick.label1.set_fontweight('bold')
ax.set_xlim([0,400])
ax.set_ylim([0,0.8])
grid(True)
show()

Красные квадраты — новые пики, обнаруженные в зеленом сигнале: введите здесь описание изображения

person Elad Joseph    schedule 09.06.2016
comment
К сожалению, я не могу определить фиксированный порог, новый обнаруженный пик может быть больше/меньше окружающих его пиков. - person Mehdi; 09.06.2016
comment
Кажется, у вас есть четко выраженные пики, которые намного сильнее шума. Просто установите пороговый уровень выше вашего уровня шума. - person Elad Joseph; 09.06.2016