0 / 0 / 1
Регистрация: 14.11.2012
Сообщений: 228
1

БПФ синусоиды

11.03.2014, 13:07. Показов 3880. Ответов 12
Метки нет (Все метки)

Делаю БПФ для массива данных (индекс-отсчет, значение массива- значение файла в момент этого отсчета), переделав его под офлайн обработку данных

Проделав это преобразование для 1 кГц синусоиды вижу бред- 2 всплеска не понятно при какой частоте. Не подскажите, пожалуйста, почему такое может быть?
Delphi
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
procedure MakeFFT;  {Процедура самого преобразования}
 var
  fftb:          TFFTBase; //класс, который реализует БПФ
  fFFTComplBuf:  ^TComplexArray;  //Буфер для хранения комплексных величин
  i: integer;
 begin
  GetMem(fFFTComplBuf, BufSize*SizeOf(TComplex)); //Выделение памяти под массив
  for i:=0 to BufSize-1 do //Заполняем данными массив
   begin
    fFFTComplBuf[i].Re := fDataBuf[i];
    fFFTComplBuf[i].Im := 0;
     end;
 
  fftb:=TFFTBase.Create(nil);
 
//FFT - выполнение БПФ
//Параметры задаваемые процедуре :
 //указатель на массив данных (с комплексными числами)
 //N - количество данных (размерность массива)
 //2^X=N (степень числа два)
 //False – прямое преобразование; True – обратное
 //Тип окна:
  // 0-прямоугольное
  // 1-треугольное
  // 2-Хэминга
  // 3-Ханна
  // 4-Блэкмана
 
  fftb.FFT(Pointer(fFFTComplBuf), BufSize, 8, False, 0);
 
  for i:=0 to BufSize-1 do //Переносим результат БПФ в исходный массив
   begin
    fDataBuf[i] := Round(fFFTComplBuf[i].Re / 500);  //заполняем массив выходными значениями (предварительно масштабируем их)
     end;
 
  fftb.Free;
  FreeMem(fFFTComplBuf, BufSize*SizeOf(TComplex)); //Освобождение памяти выделенной под массив
end;
 
procedure TForm1.Button6Click(Sender: TObject);
var i: integer;
begin
MakeFFT();   
  for i:=0 to Trunc((d/2)-1) do   
   begin
   // tmpY := Round(Buf[i]);
     Series2.addxy(i, Buf[i],'',clblue);
end;
 end;
Миниатюры
БПФ синусоиды  
__________________
Помощь в написании контрольных, курсовых и дипломных работ здесь
0
Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
11.03.2014, 13:07
Ответы с готовыми решениями:

С# БПФ
Здравствуйте, может есть у кого реализация БПФ на С# для любого количества точек?

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

БПФ и разрядность
Помогите. Нужно вичеслить розрядность (N) БП, если время преобразования равен 20нс, а частота...

Размазывание при БПФ?
Приветствую всех. Друзья, снова заезженное БПФ, с которым я совсем недавно связался и до сих пор...

12
Пишу на Delphi...иногда
1419 / 1276 / 286
Регистрация: 03.12.2012
Сообщений: 3,914
Записей в блоге: 5
11.03.2014, 18:39 2
Цитата Сообщение от Turgenev Посмотреть сообщение
вижу бред- 2 всплеска не понятно при какой частоте
все верно, так и должно быть
Простыми словами о преобразовании Фурье
спектральный анализ
1
0 / 0 / 1
Регистрация: 14.11.2012
Сообщений: 228
13.03.2014, 20:31  [ТС] 3
Цитата Сообщение от cotseec Посмотреть сообщение
все верно, так и должно быть
Простыми словами о преобразовании Фурье
спектральный анализ
Спасибо, я ознакомился с этим материалом и понял почему две дельта функции на графике (2я недействительна и при визуализации её отсекают), но все же д-функция выскакивает не на той частоте- я пересчитал относительно построенной осциллограммы частоту, которая должна быть (период синусоиды 20 отсчетов, умножаем их на величину обратную частоте дискретизации (45*10^-6 с), потом единицу делим на получившееся знаение и получаем частоту сигнала 1000 Гц). В программе же при взятии размера массива для БПФ длинной 65536 (2 в 16 степени) всплеск на 1400 Гц, при длинне массива 32768 (2 в 15 степени) всплеск на 800 Гц. Подскажите, пожалуйста, как добиться истинной частоты сигнала?
0
0 / 0 / 1
Регистрация: 14.11.2012
Сообщений: 228
16.03.2014, 18:08  [ТС] 4
никто не подскажет в связи с чем всплеск так сильно колеблется при изменении размера буфера (в коде размер задается константой BufSize), но ни разу не попадает на истинное значение 1000 Гц?
0
10218 / 6598 / 495
Регистрация: 28.12.2010
Сообщений: 21,161
Записей в блоге: 1
16.03.2014, 18:13 5
...шаг по частотам при БПФ = частота дискретизации / количество точек преобразования. Для точного определения частоты есть другие алгоритмы, к примеру Герцеля.
1
0 / 0 / 1
Регистрация: 14.11.2012
Сообщений: 228
16.03.2014, 19:12  [ТС] 6
Цитата Сообщение от raxp Посмотреть сообщение
...шаг по частотам при БПФ = частота дискретизации / количество точек преобразования.
Теперь буду знать, большое спасибо! Но я похоже недопонял, потому что при проверке (при размере буфера 65536 и частоте дискретизации 22050 Гц шаг получился 0,3, а при длине буфера 32768 шаг 0,6) шаг получился довольно мелкий все же (не сотни Герц)
[quote="raxp;5907973"Для точного определения частоты есть другие алгоритмы, к примеру Герцеля.[/quote]
За это еще одно большое спасибо! Пошел искать. Если вы знаете где можно найти пример такого преобразования на делфи- был бы вам очень признателен, если бы вы поделились ссылочкой

[size="1"]Добавлено через 10 минут[/size]
Цитата Сообщение от raxp Посмотреть сообщение
Закрепленная тема Литература по ЦОС и алгоритмам в данной ветке по ЦОС
с этим разделом ознакомился, там нужного примера более точного определения частоты звука wav файла на делфи нет((
0
10218 / 6598 / 495
Регистрация: 28.12.2010
Сообщений: 21,161
Записей в блоге: 1
16.03.2014, 19:32 7
Но я похоже недопонял, потому что при проверке (при размере буфера 65536 и частоте дискретизации 22050 Гц шаг получился 0,3, а при длине буфера 32768 шаг 0,6) шаг получился довольно мелкий все же
не размер исходного буфера, а точек преобразования. И да, шаг и будет разный.

Герцеля.
... где можно найти пример такого преобразования на делфи
тема обсосана тут Как узнать частоту звука?

А пример ДПФ по алгоритму Герцеля есть тут http://www.dsplib.ru/content/g... rtzel.html и подробно рассмотрен на Delphi, к примеру в публикации: "В.Степанов. Захват звука средствами WinAPI. - Радиолюбитель, Минск, 2014, №2, стр.12-15".

p.s.: ресурсы и исходники к нашим публикациям мы выкладываем к каждому выпуску в открытом виде http://radioliga.com/insert_2014.htm
1
2013 / 1285 / 60
Регистрация: 05.06.2010
Сообщений: 2,213
16.03.2014, 19:42 8
Алгоритм Герцеля позволяет находить один отсчет фурье, а не все сразу, как в бпф. Но он не не относится к быстрым алгоритмам, так как если считать по нему весь спектр, количество операций пропорционально 3-й степени, как в прямой формуле. Есть код на си, могу выложить если надо, не думаю что его трудно перевести на делфи или любой другой язык

Добавлено через 3 минуты
а, уже есть. И тоже на сях)

Добавлено через 3 минуты
C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
complex_t goertzel(internalType *s, int length, int idx)
{
    complex_t z;
    internalType arg = 2*pi*idx/length;
    internalType wr, wi, alpha;
    internalType v1, v2, v3;
    int i;
 
    arm_sin_cos(arg, &wi, &wr);
    alpha = 2 * wr;
    v1 = s[0];
    v2 = s[1] + alpha*v1;
    for(i=2; i<length-1; i++)
    {
        v3 = s[i] + alpha*v2 - v1;
        v1 = v2;
        v2 = v3;
    }
    v3 = s[length-1] + alpha*v2 - v1;
    z.r = v3*wr - v2,
    z.i = v3*wi;
 
    return z;
}
мой чуток оптимальнее, так как не нуждается в памяти для хранения промежуточных данных. Но в общем одно и то же. Да и что там можно выдумать нового, алгоритм то прост как молоток
1
0 / 0 / 1
Регистрация: 14.11.2012
Сообщений: 228
16.03.2014, 22:41  [ТС] 9
спасибо за ответы. Пошел по пути преобразования алгоритма на си указанного тут: http://www.dsplib.ru/content/g... rtzel.html в код на делфи. В коде по той ссылке заменил массив с данными s на свой массив с данными buf, который я получаю из wav файла (массив со значениями амплитуд в момент отсчета), а массив с промежуточными результатами v я реализую переменной длинны (но всегда равной длине массива с данными) и заполняю его нулями (в моей программе это массив fur):
Delphi
1
2
3
4
 SetLength(fur, l);
     randomize;
for c:=1 to l do
fur[c]:=random(0);
Затормозил я в конце. Не понимаю как реализовать в делфи расчет реальной и мнимой части спектрального отсчета S(1) согласно (13)(S(1) и (13) это из ссылки цитата), так, чтобы в конце концов я все же мог вывести в chart эту дельта функцию на значении частоты сигнала:
C
1
2
double SR = v[N-1]*wr-v[N-2];
    double SI = v[N-1]*wi;
Я давно общался с си и что-то вообще в ступор впал с этими выражениями. Подскажите, пожалуйста, как все же эти строчки реализовать на делфи и вывести их в chart???
Вот как я переделал из си в делфи до этого момента:
Delphi
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
procedure MakeFFT;  {Ïðîöåäóðà ñàìîãî ïðåîáðàçîâàíèÿ}
 var
 SR, SI, wi,wr ,a: double;
  N, k, i: integer;
 begin
  N:=l; //òî÷åê ÄÏÔ
     k:= 1; //íîìåð ñïåêòðàëüíîãî îòñ÷òà
 
    //ïàðàìåòð àëüôà
 a:= 2*cos(2*M_PI*k/N);
 
    //ïîâîðîòíûé ê-ò W_N^-k
     wr:=    cos(2*M_PI*k/N); //ðåàëüíàÿ ÷àñòü
     wi:=    sin(2*M_PI*k/N); //ìíèìàÿ ÷àñòü
 
    //ó÷åò v[-1] = v[-2] = 0
    fur[0]:= buf[0];
    fur[1]:= buf[1]+a*fur[0];
 
 
    i:=2;
repeat
fur[i]:= buf[i] + a*fur[i-1] - fur[i-2];
i:=i+1
until i<N;
 
    //ðåàëüíàÿ è ìíèìàÿ ÷àñòè ñïåêòðàëüíîãî îòñ÷åòà S(1) ñîãëàñíî (13)
 
 
     SR:= fur[N-1]*wr-fur[N-2];
   SI:= fur[N-1]*wi;
 
 
end;
Добавлено через 1 час 49 минут
как я думаю- нужно поместить SR и SI в цикл, в котором
будут перебираться значения массива fur и подставляться в выражения для SR и SI. Результат помещать в новый массив и в идеале в этом массиве все значения частот, кроме частоты сигнала будут иметь нулевую амплитуду (то есть квадратный корень из суммы квадратов реальной и мнимой части). Так ли это и реализуемо ли это в делфи?
0
10218 / 6598 / 495
Регистрация: 28.12.2010
Сообщений: 21,161
Записей в блоге: 1
17.03.2014, 08:46 10
http://radioliga.com/insert_20... 0sound.zip
Delphi
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
//Расчет по Герцелю
//Предварительный расчет альфа и поворотных коэффициентов для всех 882 отсчетов
 Frq:=-100;         //Начальная частота равна -100 Гц.
 For i:=0 To N-1 do begin
  k:=Trunc(N*Frq/Fs);             //Определим номер расчетной гармоники
  alpha[i]:=2*cos(2*Pi*(k/N));    //Расчет коэфф. Альфа
  wr[i] := cos(2*Pi*(k/N));       //Поворотный коэфф. реальная часть
  wi[i] := sin(2*Pi*(k/N));       //Поворотный коэфф. мнимая часть
  Frq:=Frq+F;                     //Перейдем на следующую частоту
 end;
//Цикл расчета по частотам от 50 до 10000 Гц с разносом 50 Гц
  For j:=0 to Ns-1 do begin
//Для начала расчета примем B[-1] = B[-2] = 0
   B[0] := BufData[0];
   B[1] := BufData[1]+alpha[j+1]*B[0];
//Цикл расчета B
    For i:=2 to N-1 do begin
     B[i] := BufData[i]+alpha[j+2]*B[i-1]-B[i-2];
    end;
//реальная и мнимая части спектрального отсчета
   Wr[j]:= B[N-1]*wr[j]-B[N-2];
   Wi[j]:= B[N-1]*wi[j];
//Избавимся от минусов
   Re:=Wr[j]*Wr[j];
   Im:=Wi[j]*Wi[j];
//Выделим и отмасштабируем амплитуду гармоники
   Ampl[j]:=(SQRT(Re+Im))/150000;
   If j=0 Then Ampl[j]:=0;
  end;
вот и весь алгоритм на Delphi для всего спектра.

БПФ синусоиды
1
0 / 0 / 1
Регистрация: 14.11.2012
Сообщений: 228
18.03.2014, 00:37  [ТС] 11
Цитата Сообщение от raxp Посмотреть сообщение
http://radioliga.com/insert_2014/capture%20sound.zip
Отличная ссылка! Переработал под свою прогу код по этой ссылке(просто заменив массив BufData на свой buf, в которой загружаю значения при чтении wav файла)- получил дельта-функцию на частоте всегда меньше нужной в 2 раза. Изменял все что мог в пределах разумного- нужного результата не получил. Не подскажите, пожалуйста, как и с чем я снова тут закосячил? Осталось подозрение что тип данных BufD = array[0..N-1] of Double; должен быть длинной равной длине массива buf (то есть быть переменной, в зависимости от файла). Вот мой код:
Delphi
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
 BitSample = 16;
 
  Fs = 22050;
  Fmax = 10000;           //Ìàêñèìàëüíàÿ ÷àñòîòà ñïåêòðà
  F = 50;                 //×àñòîòà ðàçíîñà ñïåêòðà
  Ns = Fmax div F;        //Êîëè÷åñòâî òî÷åê ñïåêòðà
  N = Fs div F * 2;       //Êîëè÷åñòâî îòñ÷åòîâ (ðàçìåð áóôåðà 882)
  Ch = 1;
type
//Îáúÿâèì ìàññèâû äëÿ äàííûõ
 
  BufD = array[0..N-1] of Double;
implementation
...
procedure TForm1.Button6Click(Sender: TObject);
 var
 i,j,k: integer;
  Frq: Integer;                 //Ðàñ÷èòûâàåìàÿ ÷àñòîòà
 
  B, alpha, Wr, Wi, Ampl: BufD;
  Re, Im: Double;
 begin
 Frq:=-100;         //Íà÷àëüíàÿ ÷àñòîòà ðàâíà -100 Ãö.
 For i:=0 To N-1 do begin
  k:=Trunc(N*Frq/Fs);             //Îïðåäåëèì íîìåð ðàñ÷åòíîé ãàðìîíèêè
  alpha[i]:=2*cos(2*Pi*(k/N));    //Ðàñ÷åò êîýôô. Àëüôà
  wr[i] := cos(2*Pi*(k/N));       //Ïîâîðîòíûé êîýôô. ðåàëüíàÿ ÷àñòü
  wi[i] := sin(2*Pi*(k/N));       //Ïîâîðîòíûé êîýôô. ìíèìàÿ ÷àñòü
  Frq:=Frq+F;                     //Ïåðåéäåì íà ñëåäóþùóþ ÷àñòîòó
 end;
//Öèêë ðàñ÷åòà ïî ÷àñòîòàì îò 50 äî 10000 Ãö ñ ðàçíîñîì 50 Ãö
  For j:=0 to Ns-1 do begin
//Äëÿ íà÷àëà ðàñ÷åòà ïðèìåì B[-1] = B[-2] = 0
   B[0] := buf[0];
   B[1] := buf[1]+alpha[j+1]*B[0];
//Öèêë ðàñ÷åòà B
    For i:=2 to N-1 do begin
     B[i] := buf[i]+alpha[j+2]*B[i-1]-B[i-2];
    end;
//ðåàëüíàÿ è ìíèìàÿ ÷àñòè ñïåêòðàëüíîãî îòñ÷åòà
   Wr[j]:= B[N-1]*wr[j]-B[N-2];
   Wi[j]:= B[N-1]*wi[j];
//Èçáàâèìñÿ îò ìèíóñîâ
   Re:=Wr[j]*Wr[j];
   Im:=Wi[j]*Wi[j];
//Âûäåëèì è îòìàñøòàáèðóåì àìïëèòóäó ãàðìîíèêè
   Ampl[j]:=(SQRT(Re+Im))/150000;
   If j=0 Then Ampl[j]:=0;
  end;
   For i:=0 to Ns-1 do begin
//  If Ampl[i]<2 Then Ampl[i]:=0.5;
  Chart2.Series[0].AddXY(i*50,Ampl[i]);
end;
end;
0
10218 / 6598 / 495
Регистрация: 28.12.2010
Сообщений: 21,161
Записей в блоге: 1
18.03.2014, 08:44 12
Cсылка говорите отличная, а еще ранее постом (#7) ее плохо было видно? Да, у вас есть ошибка при копипасте, посмотрите пример еще раз.
0
0 / 0 / 1
Регистрация: 14.11.2012
Сообщений: 228
18.03.2014, 11:21  [ТС] 13
Цитата Сообщение от raxp Посмотреть сообщение
Cсылка говорите отличная, а еще ранее постом (#7) ее плохо было видно?
Мой косяк, извините
Цитата Сообщение от raxp Посмотреть сообщение
Да, у вас есть ошибка при копипасте, посмотрите пример еще раз.
Я друг напротив друга вывел 2 окна- первое окно моя прога, второе- ваш пример, что сюда выкидывали и потом то, что было в примере по ссылке. Все один в один за исключением массива данных BufData (у меня он просто buf). Но после того как не нашел ошибки при копиравании побаловался с выводом частоты на график
Delphi
1
Chart2.Series[0].AddXY(i*50,Ampl[i]);
заменив i*50 на 100 все получилось, но что то мне подсказывает что так дела тут не делаются. Если вам не сложно, не могли бы вы, пожалуйста, все же указать на ошибку?
0
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
18.03.2014, 11:21

Граф-бабочка БПФ
Добрый день, в общем, такой вопрос. Дан алгоритм БПФ с прореживанием по частоте, необходимо...

БПФ на процессоре обработки сигналов
Здравствуйте! Интересует вопросы о реализации БПФ на процессоре обработки сигналов таких как TMS...

БПФ для датчика звука
Здравствуйте! Мне нужно: 1)Датчик звука, микрофон 2)Я подключаю его к компу 3)Звук, приходящий...

Правильно ли построен график БПФ?
Ребят, вопрос состоит в следующем: мне нужно понять, правильный ли у меня получился график. Что я...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
13
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.