Наши проекты:
Журнал · Discuz!ML · Wiki · DRKB · Помощь проекту |
||
ПРАВИЛА | FAQ | Помощь | Поиск | Участники | Календарь | Избранное | RSS |
[18.118.9.146] |
|
Сообщ.
#1
,
|
|
|
Доброго времени суток. Разбираюсь с FFT и IFFT. Уже какой день парюсь над обратным преобразованием Фурье. Добрые люди, может объясните в чем может быть проблема?
Звуковой сигнал 2 кГц, который имитирует программка справа, поступает на микрофон. Данные заносятся в массив и рисуется входной сигнал - верхняя диаграмма. Этот массив с данными раскладываю в FFT, получаю массив fft и рисую диаграмму - второй график сверху. Проверял на разных частотах, все здесь правильно. И наконец, я делаю обратное FFT, взяв массив fft. В итоге, 3-й график, вроде и частота сигнала та же, что и на входе, но почему-то блин амплитуда сигнала падает. И видно,что на стыке двух спектров (основного и зеркального) амплитуда минимальна. Делал эксперимент, если в спектре четко только одна гармоника, то при IFFT все нормально, амплитуда сигнала равномерна по всей длине, как на входном графике. Как только в спектре есть две гармоники и более, то сразу же график как на картинке. Как будто перемножение что ли где-то, или что может быть? Буду благодарен знатокам! procedure TFFTBase.InitMem(WindowSize: Integer); begin fWindowSize:=WindowSize; tTlbSize:=(2*fWindowSize+1)*SizeOf(Single); fSinTbl:=AllocMem(tTlbSize); fCosTbl:=AllocMem(tTlbSize); end; procedure TFFTBase.InitSinCosTbl; var i :integer; begin for i :=1 to 2*fWindowSize do begin fSinTbl[i] := (-1)*Sin(PI/i); fCosTbl[i] := Cos(PI/i); end; end; NV2 :=N shr 1; NM1 :=N-1; J :=1; if Inverse then for i :=0 to N-1 do TComplexArray(F^)[i].Im :=-TComplexArray(F^)[i].Im; for I :=1 to NM1 do begin if I<J then begin T :=TComplexArray(F^)[J-1]; TComplexArray(F^)[J-1] :=TComplexArray(F^)[I-1]; TComplexArray(F^)[I-1] :=T; end; K :=NV2; while K < J do begin J :=J - K; K :=K shr 1; end; J :=J + K; end; for L :=1 to M do begin LE :=2 shl (L-1); LE1 :=LE shr 1; U.Re :=1.0; U.Im :=0.0; W.Re :=fCosTbl[LE1]; W.Im :=fSinTbl[LE1]; for J :=1 to LE1 do begin I :=J; while I <= N do begin IP :=I + LE1; T.Re :=TComplexArray(F^)[IP-1].Re * U.Re - TComplexArray(F^)[IP-1].Im * U.Im; T.Im :=TComplexArray(F^)[IP-1].Re * U.Im + TComplexArray(F^)[IP-1].Im * U.Re; TComplexArray(F^)[IP-1].Re :=TComplexArray(F^)[I-1].Re - T.Re; TComplexArray(F^)[IP-1].Im :=TComplexArray(F^)[I-1].Im - T.Im; TComplexArray(F^)[I-1].Re:=TComplexArray(F^)[I-1].Re+T.Re; TComplexArray(F^)[I-1].Im:=TComplexArray(F^)[I-1].Im+T.Im; Inc(I,LE); end; Uo :=U; U.Re :=(Uo.Re * W.Re) - (Uo.Im * W.Im); U.Im :=(Uo.Re * W.Im) + (Uo.Im * W.Re); end; end; ImDummy :=1/Sqrt(N); if Inverse then ImDummy :=-ImDummy; ReDummy :=Abs(ImDummy); for I :=1 to N do begin TComplexArray(F^)[i-1].Re :=TComplexArray(F^)[i-1].Re*ReDummy; TComplexArray(F^)[i-1].Im :=TComplexArray(F^)[i-1].Im*ImDummy; end; DelMem; end; Прикреплённый файл111.png (53,48 Кбайт, скачиваний: 557) |
Сообщ.
#2
,
|
|
|
Здесь не приведен заголовок процедуры FFT. Как она применяется? Каков размер массива?
Какие частоты вызывают такой эффект? Получается ли то же самое с достоверно рабочим кодом FFT? |
Сообщ.
#3
,
|
|
|
Цитата MBo @ Применен типовой код: unit FFTBase; interface uses Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs, FFTFilterConst; {$R-} type TFFTBase = class(TComponent) private { Private declarations } fWindowSize, tTlbSize: Integer; fSinTbl, fCosTbl : PTbl; procedure InitSinCosTbl; protected { Protected declarations } public { Public declarations } procedure FFT( var F :Pointer; N, M :integer; const Inverse :boolean; const Window :Integer ); //N -размер F, M - N =2^M procedure simpleFFT( var F :Pointer; N, M :integer; const Inverse :boolean; const Window :Integer ); //N -размер F, M - N =2^M procedure InitMem(WindowSize: Integer); procedure DelMem; published { Published declarations } end; procedure Register; implementation {+private++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++} procedure TFFTBase.InitSinCosTbl; var i :integer; begin for i :=1 to 2*fWindowSize do begin fSinTbl[i] := (-1)*Sin(PI/i); fCosTbl[i] := Cos(PI/i); end; end; {-private----------------------------------------------------------------------} {+PUBLIC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++} procedure TFFTBase.FFT( var F :Pointer; N, M :integer; const Inverse :boolean; const Window :Integer ); //N -размер F, M - N =2^M var U, Uo, W, T :TComplex; I,IP,J,K,L,LE,LE1,NV2,NM1 :integer; ReDummy, ImDummy :extended; procedure SetWindow(const Window: Integer); var i: Integer; begin if Not Inverse then //Если не обратный case Window of //0-прямоугольное//1-треугольное//2-Хэминга//3-Ханна//4-Блэкмана 0: begin //0-прямоугольное {for i :=0 to N-1 do F[i].Re :=F[i].Re*1.2} end; 1: begin //1-треугольное for i :=0 to N-1 do begin if (i <= N div 2) then TComplexArray(F^)[i].Re :=TComplexArray(F^)[i].Re*i/(N/2) else TComplexArray(F^)[i].Re :=TComplexArray(F^)[i].Re*(N-i)/(N/2); end; end; 2: begin //2-Хэминга for i :=0 to N-1 do TComplexArray(F^)[i].Re :=TComplexArray(F^)[i].Re*(0.54-0.46*cos(2*Pi*i/N)) end; 3: begin //3-Ханна for i :=0 to N-1 do TComplexArray(F^)[i].Re :=TComplexArray(F^)[i].Re*(1+cos(2*Pi*i/N))/2; end; 4: begin //4-Блэкмана for i :=0 to N-1 do TComplexArray(F^)[i].Re :=TComplexArray(F^)[i].Re*(0.42+0.5*cos(2*Pi*i/N)+0.08*cos(4*Pi*i/N)) end; end; end; begin InitMem(N); InitSinCosTbl; SetWindow(Window); NV2 :=N shr 1; NM1 :=N-1; J :=1; if Inverse then for i :=0 to N-1 do TComplexArray(F^)[i].Im :=-TComplexArray(F^)[i].Im; for I :=1 to NM1 do begin if I<J then begin T :=TComplexArray(F^)[J-1]; TComplexArray(F^)[J-1] :=TComplexArray(F^)[I-1]; TComplexArray(F^)[I-1] :=T; end; K :=NV2; while K < J do begin J :=J - K; K :=K shr 1; end; J :=J + K; end; for L :=1 to M do begin LE :=2 shl (L-1); LE1 :=LE shr 1; U.Re :=1.0; U.Im :=0.0; W.Re :=fCosTbl[LE1]; W.Im :=fSinTbl[LE1]; for J :=1 to LE1 do begin I :=J; while I <= N do begin IP :=I + LE1; T.Re :=TComplexArray(F^)[IP-1].Re * U.Re - TComplexArray(F^)[IP-1].Im * U.Im; T.Im :=TComplexArray(F^)[IP-1].Re * U.Im + TComplexArray(F^)[IP-1].Im * U.Re; TComplexArray(F^)[IP-1].Re :=TComplexArray(F^)[I-1].Re - T.Re; TComplexArray(F^)[IP-1].Im :=TComplexArray(F^)[I-1].Im - T.Im; TComplexArray(F^)[I-1].Re:=TComplexArray(F^)[I-1].Re+T.Re; TComplexArray(F^)[I-1].Im:=TComplexArray(F^)[I-1].Im+T.Im; Inc(I,LE); end; Uo :=U; U.Re :=(Uo.Re * W.Re) - (Uo.Im * W.Im); U.Im :=(Uo.Re * W.Im) + (Uo.Im * W.Re); end; end; ImDummy :=1/Sqrt(N); if Inverse then ImDummy :=-ImDummy; ReDummy :=Abs(ImDummy); for I :=1 to N do begin TComplexArray(F^)[i-1].Re :=TComplexArray(F^)[i-1].Re*ReDummy; TComplexArray(F^)[i-1].Im :=TComplexArray(F^)[i-1].Im*ImDummy; end; if not Inverse then //Вычисление Мощности С=SQRT(А^2+B^2) при -прямом- БПФ for i :=0 to N-1 do TComplexArray(F^)[i].Re :=sqrt(sqr(TComplexArray(F^)[i].Re) +sqr(TComplexArray(F^)[i].Im)); DelMem; //Освобождение ресурсов end; procedure TFFTBase.InitMem(WindowSize: Integer); begin fWindowSize:=WindowSize; tTlbSize:=(2*fWindowSize+1)*SizeOf(Single); fSinTbl:=AllocMem(tTlbSize); fCosTbl:=AllocMem(tTlbSize); end; procedure TFFTBase.DelMem; begin FreeMem(fCosTbl, tTlbSize); FreeMem(fSinTbl, tTlbSize); end; Так происходит при любых частотах. BufSize = N = 1024. Другие FFT еще не пробовал... Из моей картинки, как я понял, на втором и третьем графиках только половины сигналов спектра и IFFT, правые части - зеркальны. Восстановить целиком сигнал не выйдет, получается?... Запускается этот FFTBase отсюда, после комментария // 4-Блэкмана: procedure MakeIFFT; //Делаем IFFT 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, 10, True, 2); for i:=0 to BufSize-1 do //Переносим результат БПФ в исходный массив begin //заполняем массив выходными значениями (предварительно масштабируем их) fDataBuf[i] := Round(fFFTComplBuf[i].Re / 1); end; fftb.Free; FreeMem(fFFTComplBuf, BufSize*SizeOf(TComplex)); //Освобождение памяти выделенной под массив end; И выполняются процедуры в первом коде по очереди: InitMem(N); InitSinCosTbl; SetWindow(Window); затем само Фурье - NV2 :=N shr 1; NM1 :=N-1; J :=1; и т.д. Результат переносится в fDataBuf[i] во втором коде. |
Сообщ.
#4
,
|
|
|
Так же на 10 кГц, видно на картинке...
На всех частотах. Однако, к примеру на 5 кГц, я в массиве обнулил небольшие гармоники, оставил только одну - самую большую - все нормально, сигнал на выходе такой же как на входе. Может что из-за окна? При FFT применяю окно 2-Хэмминга. Очень похоже. А при IFFT окна не применяются, исходя из кода. Прикреплённый файл222.png (40,23 Кбайт, скачиваний: 639) Прикреплённый файл333.png (49,93 Кбайт, скачиваний: 620) |
Сообщ.
#5
,
|
|
|
Сделал сигнал длиной 1024 отсчета, ввел одну, две, три частоты.
Не наблюдаю такого эффекта как с окном, так и без окна. Может быть, выводятся графики не того, что есть на самом деле? var Re, Im: array of Single; i, N: Integer; Amp: Double; begin N := 1024; SetLength(Re, N); SetLength(Im, N); for i := 0 to N - 1 do begin Re[i] := Sin(50 * i * 2 * Pi / N) + 0.3 * Sin(67 * i * 2 * Pi / N) + 0.3 * Sin(91 * i * 2 * Pi / N); //окно, его можно отключить Re[i] := Re[i] * (0.54-0.46*cos(2*Pi*i/N)); Im[i] := 0; Series1.AddXY(i, Re[i]); end; FFT1D(Re, Im, N, 1); for I := 0 to N - 1 do begin Amp := Hypot(Re[i], Im[i]); Series2.AddXY(i, Amp); end; FFT1D(Re, Im, N, 0); //обратное for I := 0 to N - 1 do begin Series4.AddXY(i, Re[i] + 3.1); //сдвинем для удобства end; Прикреплённый файлfft.png (13,99 Кбайт, скачиваний: 597) |
Сообщ.
#6
,
|
|
|
Скорее всего у меня выводится сигнал зеркально и выходит, что левая и правая половинки вместе образуют по середине провал по амплитуде... У Вас амплитуда по краям минимальная, а у меня наоборот - посередине.
И все же, если поставить окно прямоугольного типа, то эти провалы куда меньше. Найду другой FFT проверю как там. Спасибо. |
Сообщ.
#7
,
|
|
|
Амплитуда минимальная по краям из-за применения окна
|