Отсечение БПФ Матрица - PullRequest
Новая
       11

Отсечение БПФ Матрица

2 голосов
/ 01 июня 2009

Обработка звука довольно новая для меня. И в настоящее время использует Python Numpy для обработки волновых файлов. После вычисления матрицы БПФ я получаю значения мощности шума для несуществующих частот. Я заинтересован в визуализации данных и точность не является высоким приоритетом. Есть ли безопасный способ вычислить значение отсечения для удаления этих значений, или я должен использовать все матрицы FFT для каждого набора выборок, чтобы получить среднее число?

привет

Edit: 

    from numpy import *
    import wave
    import pymedia.audio.sound as sound
    import time, struct
    from pylab import ion, plot, draw, show

    fp = wave.open("500-200f.wav", "rb")
    sample_rate = fp.getframerate()
    total_num_samps = fp.getnframes()
    fft_length = 2048.
    num_fft = (total_num_samps / fft_length ) - 2
    temp = zeros((num_fft,fft_length), float)

    for i in range(num_fft):
        tempb = fp.readframes(fft_length);
        data = struct.unpack("%dH"%(fft_length), tempb)
        temp[i,:] = array(data, short) 
    pts = fft_length/2+1
    data = (abs(fft.rfft(temp, fft_length)) / (pts))[:pts]

    x_axis = arange(pts)*sample_rate*.5/pts
    spec_range = pts
    plot(x_axis, data[0])
    show()

Вот график в нелогарифмическом масштабе для файла синтетических волн, содержащего синусоидальный сигнал 500 Гц (затухание) + 200 Гц, созданный с использованием Goldwave.

image

Ответы [ 4 ]

3 голосов
/ 01 июня 2009

Имитированные осциллограммы не должны отображать БПФ, подобные вашей фигуре, поэтому что-то очень неправильно, и, вероятно, не с БПФ, а с входной осциллограммой. Основная проблема в вашем сюжете - это не пульсации, а гармоники около 1000 Гц и субгармоника при 500 Гц. Имитированная форма волны не должна показывать ничего из этого (например, см. Мой график ниже).

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

Мне удалось получить довольно близкий дубликат к вашему БПФ, применив отсечение к сигналу, как было предложено как для субгармонических, так и для высших гармоник (и Тревора). Вы можете ввести обрезку либо в симуляции, либо в распаковке. В любом случае, я обошел это путем создания волновых форм в numpy для начала.

Вот как должно выглядеть правильное БПФ (т.е. в основном идеально, за исключением расширения пиков из-за окон)

alt text

Вот один из осциллограмм, которые были обрезаны (и очень похож на ваше БПФ, от субгармоники до точного паттерна трех высших гармоник около 1000 Гц)

alt text Вот код, который я использовал для генерации этих 

from numpy import *
from pylab import ion, plot, draw, show, xlabel, ylabel, figure

sample_rate = 20000.
times = arange(0, 10., 1./sample_rate)
wfm0 = sin(2*pi*200.*times)
wfm1 = sin(2*pi*500.*times) *(10.-times)/10.
wfm = wfm0+wfm1
#  int test
#wfm *= 2**8
#wfm = wfm.astype(int16)
#wfm = wfm.astype(float)
#  abs test
#wfm = abs(wfm)
#  clip test
#wfm = clip(wfm,  -1.2, 1.2)

fft_length = 5*2048.
total_num_samps = len(times)
num_fft = (total_num_samps / fft_length ) - 2
temp = zeros((num_fft,fft_length), float)

for i in range(num_fft):
    temp[i,:] = wfm[i*fft_length:(i+1)*fft_length] 
pts = fft_length/2+1
data = (abs(fft.rfft(temp, fft_length)) / (pts))[:pts]

x_axis = arange(pts)*sample_rate*.5/pts
spec_range = pts
plot(x_axis, data[2], linewidth=3)
xlabel("freq (Hz)")
ylabel('abs(FFT)')
show()
2 голосов
/ 01 июня 2009

БПФ, потому что они оконные и дискретизированные вызывают алиасинг и выборку также в частотной области. Фильтрация во временной области - это просто умножение в частотной области, поэтому вы можете просто применить фильтр, который просто умножает каждую частоту на значение для функции для фильтра, который вы используете. Например, умножьте на 1 в полосе пропускания и на ноль все остальные. Неожиданные значения, вероятно, вызваны псевдонимами, где более высокие частоты сворачиваются до тех, которые вы видите. Исходный сигнал должен быть ограничен полосой до половины вашей частоты дискретизации , иначе вы получите aliasing . Еще большее беспокойство вызывает алиасинг, который искажает интересующую область, потому что для этой полосы частот вы хотите знать, что частота соответствует ожидаемой.

Другая вещь, которую нужно иметь в виду, заключается в том, что когда вы извлекаете часть данных из волнового файла, вы математически умножаете его на прямоугольную волну. Это приводит к свертке sinx / x с частотной характеристикой, чтобы минимизировать это, вы можете умножить исходный оконный сигнал на что-то вроде окна Ханнинга .

1 голос
/ 01 июня 2009

Я не знаю достаточно из вашего вопроса, чтобы на самом деле ответить на что-то конкретное.

Но вот несколько вещей из моего собственного опыта написания БПФ:

  • Убедитесь, что вы следуете правилу Найквиста
  • Если вы просматриваете линейный выход БПФ ... у вас будут проблемы с отображением собственного сигнала и вы думаете, что все сломано. Убедитесь, что вы смотрите на дБ вашей величины БПФ. (то есть "сюжет (10 * log10 (abs (fft (x))))")
  • Создайте unitTest для вашей функции FFT (), подавая сгенерированные данные как чистый тон. Затем передайте те же самые сгенерированные данные в FFT () Matlab. Сделайте разность абсолютных значений между двумя сериями выходных данных и убедитесь, что максимальная разница абсолютных значений составляет что-то вроде 10 ^ -6 (т.е. единственная разница вызвана небольшими ошибками с плавающей запятой)
  • Убедитесь, что вы открываете свои данные

Если все эти три вещи сработают, то у вас все в порядке. И ваши входные данные, вероятно, проблема.

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

1 голос
/ 01 июня 2009

Стоит отметить, что для 1D БПФ первый элемент (индекс [0]) содержит член постоянного тока (нулевой частоты), элементы [1:N/2] содержат положительные частоты, а элементы [N/2+1:N-1] содержат отрицательные частоты. Поскольку вы не предоставили пример кода или дополнительную информацию о выходе вашего БПФ, я не могу исключить возможность того, что «значения мощности шума на несуществующих частотах» - это не просто отрицательные частоты вашего спектра.

РЕДАКТИРОВАТЬ : Здесь - это пример БПФ radix-2, реализованного на чистом Python, с простой процедурой тестирования, которая находит БПФ прямоугольного импульса [1.,1.,1.,1.,0.,0.,0.,0.]. Вы можете запустить пример на codepad и увидеть, что БПФ этой последовательности 

[0j,                    Negative frequencies
(1+0.414213562373j),    ^
0j,                     |
(1+2.41421356237j),     |
(4+0j),                <= DC term
(1-2.41421356237j),     |
0j,                     v
(1-0.414213562373j)]    Positive frequencies

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

...