191 / 127 / 52
Регистрация: 19.01.2010
Сообщений: 518
1

Matlab vs FFTW. FFT

22.11.2013, 00:20. Показов 4033. Ответов 16
Метки нет (Все метки)

Приветствую. Делаю программу, отображающую спектр wav-файла. Столкнулся с кое-каким непониманием: при вычислении БПФ в матлабе получаются одни значения, а в своей программе (использую библиотеку fftw) - другие. Проверял правильность работы функции из этой библиотеки, считая ДПФ небольшой последовательности в маткаде - все верно.
Значения амплитуд из wav-файла считываю правильно, сверялся с тем же матлабом (функция wavread). Может надо как-то нормализовывать результат БПФ? Гуглил, читал, что каждый коэффициент БПФ делят на N, но результат все равно не совпадает с матлабовским.
Может кто знает почему так получается, или это несовпадение нормально?
0
Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
22.11.2013, 00:20
Ответы с готовыми решениями:

Параллелизация FFT на Matlab
Доброго времени суток! Прошу помощи по такому вопросу: Есть матрица NxM, нужно посчитать Фурье в...

Использование FFT в matlab
Здравсвтвуйте. Cтолкнулся со след проблемой.вычисляется fft в след коде log2fftlen = 9;...

Настройка блока FFT MatLab
Требуется собрать схему(см.скрин 1). Получилось собрать схему(см.скрин 2), сама схема собралась...

FFTW
Доброго времени суток! Товарищи, подскажите, работал ли кто-нибудь с бибилиотекой FFTW в Qt...

16
10218 / 6598 / 495
Регистрация: 28.12.2010
Сообщений: 21,161
Записей в блоге: 1
22.11.2013, 10:38 2
...приведите картинки для сравнения и весь сопутствующий код, догадываться где вы или не вы накосячили нет особого желания.

делят на N,
делят.
1
191 / 127 / 52
Регистрация: 19.01.2010
Сообщений: 518
22.11.2013, 13:24  [ТС] 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
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
55
56
57
58
TStringList *SL= new TStringList;
if(OpenDialog1->Execute())
        Fname = OpenDialog1->FileName;
SL->LoadFromFile(Fname);             // Загружаем файл со значениями амплитуд wav
AnsiString    type          = SL->Strings[0];      // В 1 строке файла инф. о типе (моно\стерео)                  
unsigned int SampleRate = StrToInt(SL->Strings[1]);  // Во 2 строке инф. о частоте дискретизации
unsigned int size           = SL->Count - 2;       // Количество аплитуд (N)
 
unsigned int half_size = floor(size/2) + 1;   // N/2+1
 
double *source = new double[size];
double *result = new double[size];
 
/* Загоняем в массив значения амплитуд wav-файла */
for(i = 0; i < size; i++)
        source[i] = StrToFloat(SL->Strings[i + 2]);
 
/* Выделяем память под массив для сохранения комплексного результата */
fftw_complex  *result_c = (fftw_complex*) fftw_malloc(size * sizeof(fftw_complex));
/* Построение плана вычисления БПФ. На входе: действительные. На выходе: комплексные значения */
fftw_plan     plan      = fftw_plan_dft_r2c_1d(size, source, result_c, FFTW_ESTIMATE);
/* Выполнение плана (сам расчет БПФ) */
fftw_execute(plan);
 
TStringList *out= new TStringList;
double Re2, Im2;
 
/* Считаем модуль бпф */
for(i = 0; i < half_size; i++)
{
       Re2 = result_c[i][0] * result_c[i][0];
       Im2 = result_c[i][1] * result_c[i][1];
       result[i] = sqrt(Re2 + Im2) / size;
       out->Add(FloatToStr(result[i]));     // Формируем запись о результате бпф
}
 
/* Сохраняем запись в текстовик (по этой записи я и ориентируюсь о результате вычисления БПФ) */
out->SaveToFile("out.txt");    
 
/* Формируем массив частот для оси абсцисс */
double      *freq = new double[size];
for(i = 0; i < half_size; i++)
{
      freq[i] = i + 1;
      freq[i] = SampleRate/2 * freq[i] / half_size;
}
 
delete out;
 
/* Рисуем график */
for(i = 1; i < half_size; i++)
      Series1->AddXY(freq[i], result[i], clRed);
 
fftw_free(result_c);
fftw_destroy_plan(plan);
delete []freq;
delete []source;
delete []result;

Мой график:
Кликните здесь для просмотра всего текста
Matlab vs FFTW. FFT

Код матлаб (брал отсюда) и график:
Кликните здесь для просмотра всего текста
Matlab M
1
2
3
4
5
6
[s, fs] = wavread('D:\1-mono.wav');
spec = abs(fft(s));
spec1 = spec(1:floor(end/2));
freq = (fs/2)*(1:length(spec1))/length(spec1);
plot(freq, spec1);
xlabel('Frequency: Hz'); ylabel('Amplitude spectrum');
Matlab vs FFTW. FFT
1
2013 / 1285 / 61
Регистрация: 05.06.2010
Сообщений: 2,213
22.11.2013, 14:19 4
матлаб как раз использует внутри функции fft() библиотеку fftw, так что работают они одинаково. Ищите у себя ошибку. Коэффициенты должны совпадать без всяких делений на N.
1
10218 / 6598 / 495
Регистрация: 28.12.2010
Сообщений: 21,161
Записей в блоге: 1
22.11.2013, 14:31 5
1- " result[i] = sqrt(Re2 + Im2) / size;" приведите к единообразию вашу с матлабовской осью по амплитуде, чтобы и у вас были децибеллы, возьмите дополнительно десятичный логарифм от result-а и умножьте на 20
2- что-то с отображением значений по оси 0Y в вашей программе на графике, не видно сколько реально, везде 0.001 по оси
3- что за манипуляции с i? Она и без этого инкрементируется:
C++
1
2
3
4
5
for(i = 0; i < half_size; i++)
{
      freq[i] = i + 1;
      freq[i] = SampleRate/2 * freq[i] / half_size;
}
У меня:
Delphi
1
2
3
for i := 0 to length(outwav.dy) - 1 do begin //
    freq[i]:= i * nSamplesPerSec {44100} * nChannels {2 - стерео и 1 - моно} / (2 * length(outwav.dy));
...
1
2013 / 1285 / 61
Регистрация: 05.06.2010
Сообщений: 2,213
22.11.2013, 14:39 6
raxp, пункт 2 верно, возможно это единственная проблема. В матлабе он тоже считает амплитуды без перевода в децибелы. В 3 пункте freq = i+1 ничего дополнительно не инкрементируется вы поторопились, там просто смещается на единицу. Предлагаю для начала проверить работу функции на короткой последовательности например 4 числа и сравнить результаты(числа а не графики)
1
10218 / 6598 / 495
Регистрация: 28.12.2010
Сообщений: 21,161
Записей в блоге: 1
22.11.2013, 15:05 7
без перевода в децибелы
мне удобнее смотреть на нормированный график, а код матлаба действительно не смотрел.

смещается на единицу
этот дополнительный инкремент не нужен, он уже в самом цикле задается через i++:
C++
1
for(инициализация; проверка условия; изменение) оператор;
1
2013 / 1285 / 61
Регистрация: 05.06.2010
Сообщений: 2,213
22.11.2013, 15:27 8
Цитата Сообщение от raxp Посмотреть сообщение
дополнительный инкремент не нужен
действительно не нужен, возможно Selot для совместимости с матлабовской индексацией сместил на единицу отсчеты частоты. Это не инкремент. В децибелах действительно удобнее, мне тоже так привычнее.
1
191 / 127 / 52
Регистрация: 19.01.2010
Сообщений: 518
22.11.2013, 15:28  [ТС] 9
Цитата Сообщение от vital792 Посмотреть сообщение
Предлагаю для начала проверить работу функции на короткой последовательности например 4 числа и сравнить результаты(числа а не графики)
Проверял, все совпадает.
Цитата Сообщение от vital792 Посмотреть сообщение
матлаб как раз использует внутри функции fft() библиотеку fftw, так что работают они одинаково. Ищите у себя ошибку. Коэффициенты должны совпадать без всяких делений на N.
Убрал деление на N, вот график с дБ:
Matlab vs FFTW. FFT

И обычный:
Matlab vs FFTW. FFT


freq[i] = i + 1; изменил на freq[i] = i;
0
2013 / 1285 / 61
Регистрация: 05.06.2010
Сообщений: 2,213
22.11.2013, 15:33 10
Цитата Сообщение от Selot Посмотреть сообщение
Проверял, все совпадает
так выходит проблема только в отображении? Или уже нет проблем?
0
191 / 127 / 52
Регистрация: 19.01.2010
Сообщений: 518
22.11.2013, 15:53  [ТС] 11
Цитата Сообщение от vital792 Посмотреть сообщение
так выходит проблема только в отображении? Или уже нет проблем?
Проблема осталась. Вот, смотрите...
Беру небольшую последовательность из 8 чисел: 40, 23, -2, 0, 0.2, 11, -6, 34.
Загоняю их в матлаб и считаю бпф:
Matlab vs FFTW. FFT

Загоняю их в свою прогу и считаю бпф:
Matlab vs FFTW. FFT

Все круто, считает правильно.

Теперь беру амплитуды из wav и загоняю в матлаб:
Matlab vs FFTW. FFT

А теперь в свою прогу:
Matlab vs FFTW. FFT

И даже близко не совпадает.
0
10218 / 6598 / 495
Регистрация: 28.12.2010
Сообщений: 21,161
Записей в блоге: 1
22.11.2013, 15:58 12
Приложите ваш файл со значениями амплитуд WAV, погляжу у себя на своих FFTW.
1
191 / 127 / 52
Регистрация: 19.01.2010
Сообщений: 518
22.11.2013, 16:05  [ТС] 13
Тут ограничение на размер тхт файлов стоит, поэтому пришлось в его архив:audio_data.rar
0
2013 / 1285 / 61
Регистрация: 05.06.2010
Сообщений: 2,213
22.11.2013, 16:54 14
для ваших данных код на матлабе:
Matlab M
1
2
s = load('audio_data.txt');
S=fft(s);
и код на си (без плюсов, мне так удобнее):
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
25
26
27
28
29
#include "fftw3.h"
#include <stdio.h>
#include <malloc.h>
 
int main()
{
    int i, L = 8313;
    double *signal;
    fftw_complex *spec;
    fftw_plan p;
    FILE *f = fopen("audio_data.txt", "r");
 
    signal = (double*)calloc(L, sizeof(double));
    spec = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * L);
 
    for(i=0; i<L; ++i)
    {
        fscanf(f, "%lf\n", &signal[i]);
    }
 
    p = fftw_plan_dft_r2c_1d(L, signal, spec, FFTW_ESTIMATE);
    fftw_execute(p);
 
    fftw_destroy_plan(p);
    fftw_free(spec);
    free(signal);
    fclose(f);
    return 0;
}
результат одинаков (да оно и понятно, другого и не ожидалось). Возможно ошибка у вас на этапе чтения wav. Залейте и его, попробую с ним.
1
191 / 127 / 52
Регистрация: 19.01.2010
Сообщений: 518
22.11.2013, 17:42  [ТС] 15
1-mono.rar
0
191 / 127 / 52
Регистрация: 19.01.2010
Сообщений: 518
22.11.2013, 17:56  [ТС] 16
Спасибо всем, товарищи. Теперь все работает! И графики 1 в 1 Проблема была с точностью считывании амплитуд из wav.
1
2013 / 1285 / 61
Регистрация: 05.06.2010
Сообщений: 2,213
22.11.2013, 18:25 17
Matlab M
1
2
3
4
5
6
7
8
9
clear; clc; close all;
 
s=wavread('1-mono.wav');
S = fft(s);
plot(abs(S))
SS = load('spec.txt');
hold on;
% figure;
plot(SS, 'r')
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
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
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
#include "fftw3.h"
#include <stdio.h>
#include <malloc.h>
#include <math.h>
 
typedef struct
{
    char id_riff[4];
    long len_riff;
 
    char id_chuck[4];
    char fmt[4];
    long len_chuck;
 
    short type;
    short channels;
    long freq;
    long bytes;
    short align;
    short bits;
 
    char id_data[4];
    long len_data;
}titleWave;
 
int main()
{
    int i, L = 0;
    fftw_plan p;
    titleWave title;
    FILE *f = fopen("1-mono.wav", "rb");
    FILE *fRes = fopen("spec.txt", "w");
 
    short *s;
    double *signal;
    fftw_complex *spec;
/*
    fseek(f, 0, SEEK_END);
    L = (ftell(f) - sizeof(titleWave))/ sizeof(short);
    fseek(f, 0, SEEK_SET);
*/
    fread(&title, sizeof(title), 1, f);
    L = title.len_data / sizeof(short);
 
    s = (short*)calloc(L, sizeof(short));
    signal = (double*)calloc(L, sizeof(double));
    spec = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * L/2 +1);
 
    fread(s, sizeof(short), L, f);
 
    for(i=0; i<L; ++i)
    {
        signal[i] = (double)s[i] / (1<<15); // /32768 for normalize
    }
    fclose(f);
    free(s);
 
    p = fftw_plan_dft_r2c_1d(L, signal, spec, FFTW_ESTIMATE);
    fftw_execute(p);
    fftw_destroy_plan(p);
 
    for(i=0; i<L/2; ++i)
    {
        signal[i] = sqrt(spec[i][0] * spec[i][0] + spec[i][1] * spec[i][1]);
        fprintf(fRes, "%f, ", signal[i]);
    }
 
    fftw_free(spec);
    free(signal);
    fclose(fRes);
 
    return 0;
}
не успел. Но вобщем рад что мое предположение оказалось правильным. Удачи!
1
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
22.11.2013, 18:25

Как подключить библиотеку FFTW в C#?
Возникла проблема с подключением fftw библиотеки, если добавляю dll файл через ссылки появляется...

Как скомпилировать библиотеку FFTW?
Добрый вечер. Скажите кто знает, как собрать библиотеку FFTW? (преобразование Фурье) Скачал ее, в...

Как запустить библиотеку FFTW
Довольно хорошо разобрался в вопросе. Начнем... Нашел, как другие люди делали БПФ....

Результат БПФ при использовании FFTW
Всем доброго дня! У меня есть массив данных, над которыми я выполняю БПФ, используя библиотеку...


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

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

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin
Copyright ©2000 - 2022, CyberForum.ru