Накапливается ошибка в fft

alexandr.krupnov

Выполнил ДПФ на СИ. исходные данные:
Входной сигнал - синус. Частота сигнала 1000 Гц, Амплитуда 100 В, Частота дискретизации 32 000 Гц, количество выборок = 128.
Массив действительных чисел формирую по формуле f(t) = A*sin(2*pi*f*t), t = 0,00003125
Массив мнимых чисел зануляю
По теории шаг частоты получаю в 250 Гц, то есть коэффициент (0) имеет частоту 0 Гц - это постоянная составляющая.
Шаг частоты для коэффициентов считаю по формуле Частота дискретизации/Количество выборок
коэффициент (1) имеет частоту 250 Гц, амплитуда = 0 (по теориии)
коэффициент (2) имеет частоту 500 Гц, амплитуда = 0 (по теориии)
коэффициент (3) имеет частоту 750 Гц, амплитуда = 0 (по теориии)
коэффициент (4) имеет частоту 1 000 Гц. Здесь мы должны получить амплитуду гармоники, равную половине амплитуды исходной гармоники (A' = A/2).
По зеркальному эффекту по теории должен получить также амплитуду гармоники, равную половине амплитуды исходной гармоники (A' = A/2) на коэффициенте 128-4=124.
Так вот, ДПФ работает правильно. Выполняю ОДПФ, получаю исходный сигнал.
Делаю БПФ, в итоге имею коэффициент(4) = 59,1 В (на 9,1 вольт больше), а далее через каждые 8 выборок вместо нулей некое число, порядка 1 - 8 Вольт
В вольтах указываю так как после бпф вычисляю амплитуду.
Прилагаю код, может кто посмотрит и подскажет. Спасибо

[img]/modules/file/icons/application-octet-stream.png[/img] fft.rar

folk · Сообщение **folk** » 02 апр 2014, 18:35

А если от чистого синуса что получается?

alexandr.krupnov

folk писал(а):Source of the post
А если от чистого синуса что получается?

Если Вы имеете ввиду от сигнала f(t)=sin(2*pi*f*t), f = 1000Гц, то имею следующую картину:
на 4 коэффициенте амплитуда равна 0.591
на 12 коэффициенте А = 0.069
на 20 коэффициенте А = 0.012
на 28 коэфффициенте А = 0.028
на 36 коэффициенте А = 0.036
и так далее через каждые 8 коэффициентов. На 124 коэффициенте А = 0.362
Остальные коэффициенты занулились.

Alexu007 · Сообщение **Alexu007** » 03 апр 2014, 11:27

Что такое ДПФ, БПФ, ОДПФ?

grigoriy · Сообщение **grigoriy** » 03 апр 2014, 12:57

alexandr.krupnov

Alexu007 писал(а):Source of the post
Что такое ДПФ, БПФ, ОДПФ?

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

guryev · Сообщение **guryev** » 04 апр 2014, 13:14

alexandr.krupnov писал(а):Source of the post
Прилагаю код, может кто посмотрит и подскажет. Спасибо

Может, я что-то путаю, но похоже, у вас неправильно закодировано произведение комплексных чисел в процедуре coeff_func():

Код: Выбрать все

                  c = ReX[cycle+shift_butterfly];
                  d = ImX[cycle+shift_butterfly];
                  e = Re[step*cycle_W];
                  f = Im[step*cycle_W];
                  g = c*e;
                  h = d*f;
                  j = c*f;
                  k = d*e;
                  l = g + h;
                  m = k - j;

Насколько я понял, $$l$$

- это действительная часть произведения, а $$m$$

- мнимая. Но тогда вместо

Код: Выбрать все

 l = g + h;
 m = k - j;

должно быть

Код: Выбрать все

 l = g - h;
 m = k + j;

alexandr.krupnov

guryev писал(а):Source of the post
alexandr.krupnov писал(а):Source of the post
Прилагаю код, может кто посмотрит и подскажет. Спасибо
Может, я что-то путаю, но похоже, у вас неправильно закодировано произведение комплексных чисел в процедуре coeff_func():
Код: Выбрать все
c = ReX[cycle+shift_butterfly]; d = ImX[cycle+shift_butterfly]; e = Re[step*cycle_W]; f = Im[step*cycle_W]; g = c*e; h = d*f; j = c*f; k = d*e; l = g + h; m = k - j;

Насколько я понял, - это действительная часть произведения, а - мнимая. Но тогда вместо
Код: Выбрать все
l = g + h; m = k - j;
должно быть
Код: Выбрать все
l = g - h; m = k + j;

Вы правы l - действительная часть, а m - мнимая часть произведения в FFT. Но даже если я правлю код по Вашему предложению, не меняется ровным счётом ничего.
Там ведь по формуле
a = a+b*w
b = a-b*w,
где
w = ReW - ImW = cos \alpha - j*sin\alpha
a = ReA + ImA
b = ReB + ImB. Примем ReB = ReA[cycle+shift], ImB = ImA[cycle+shift]
тогда имеем
b*w = (ReA[cycle+shift]+j*ImA[cycle+shift])*(ReW-j*ImW) = ReA[cycle+shift]*ReW - ReA[cycle+shift]*j*ImW +ReW*j*ImA[cycle+shift] - j*ImA[cycle+shift]*j*ImW =
= ( ReA[cycle+shift]*ReW - j*j*ImA[cycle+shift]*ImW) + j*( ReW*ImA[cycle+shift] - ReA[cycle+shift]*ImW ) =
( ReA[cycle+shift]*ReW - (-1)*ImA[cycle+shift]*ImW) + j*( ReW*ImA[cycle+shift] - ReA[cycle+shift]*ImW ) =
( ReA[cycle+shift]*ReW + ImA[cycle+shift]*ImW) + j*( ReW*ImA[cycle+shift] - ReA[cycle+shift]*ImW )
Тогда в левой скобке ответа мы получим действительную часть а в правой - мнимую
И подставим в формулы
a = a+b*w
b = a-b*w,
вместо b*w
Получим
a = a+b*w = (ReA+j*ImA) + ( ReA[cycle+shift]*ReW + ImA[cycle+shift]*ImW) + j*( ReW*ImA[cycle+shift] - ReA[cycle+shift]*ImW ) =
= (ReA + ( ReA[cycle+shift]*ReW + ImA[cycle+shift]*ImW) ) + j (ImA + (ReW*ImA[cycle+shift] - ReA[cycle+shift]*ImW) )

b = a-b*w = (ReA+j*ImA) - ( ReA[cycle+shift]*ReW + ImA[cycle+shift]*ImW) + j*( ReW*ImA[cycle+shift] - ReA[cycle+shift]*ImW ) =
= (ReA - ( ReA[cycle+shift]*ReW + ImA[cycle+shift]*ImW) ) + j (ImA - (ReW*ImA[cycle+shift] - ReA[cycle+shift]*ImW) )

Посмотрите пожалуйста, голова уже не видит ошибки

alexandr.krupnov

Кстати, я проследил такую закономерность:
Какую бы амплитуду я не вводил, на выходе получается не половина её как должно быть по теории, а ровно в 1,1834 раза больше.
Например А = 1, тогда после fft я должен получить 0,5. А я получаю 0,59117 (в 1,1834 раза)
ввожу А = 30, должен получить 15 но получаю 17,7522 (в 1,1834 раза) и так далее

alexandr.krupnov

И ещё одно наблюдение.
Если частота исходного сигнала 1000 Гц, то ошибка на +18% (в 1,18... раз)
если частота 2000 Гц, то ошибка на +9,6%
если 3000 Гц то -6,8%
если 4000 Гц то ошибки нет
если 5000 Гц то -4,4%
если 6000 Гц то +0,02%
если 7000 Гц то -11%
если 8000 Гц то ошибки нет
если 9000 Гц то -9%
если 10000 Гц то -3,9%
если 11000 Гц то -13,7%
если 12000 Гц то ошибки нет

yourdomain.com