FFT算法实现
FFT算法實現 Delphi / Windows SDK/API
http://www.delphi2007.net/DelphiBase/html/delphi_20061210124307209.html
type ? ? TComplexData ? = ? record ? //自定義復數類型 ?
? ? ? ? ? ? vReal,vImaginary: ? Double; ?
? end; ?
? TComplexDataArray ? = ? array ? of ? TComplexData; ? //實數數組 ? ?
? TExtendedArray ? = ? array ? of ? Extended; ? //復數數組 ?
? ?
? function ? ComplexAdd(c1,c2:TComplexData):TComplexData; ? ? ? ? ? ? ? //復數相加 ?
? function ? ComplexMinus(c1,c2:TComplexData):TComplexData; ? ? ? ? ? //復數相減 ?
? function ? ComplexMultiply(c1,c2:TComplexData):TComplexData; ? ? //復數相乘 ?
? procedure ? ComplexArrayWrap(var ? a1,a2:TComplexDataArray); ? ? ? ? //復數數組互換 ?
? procedure ? DoFFT(TD:TExtendedArray; ? var ? FD:TExtendedArray; ? r:Integer); ?
? ?
? procedure ? Tmainform.DoFFT(TD:TExtendedArray; ? var ? FD:TExtendedArray; ? r:Integer); ?
? var ?
? ? ? ? ? ? ? ? ? i,j,k ? : ? Integer; ? //循環變量 ?
? ? ? ? ? ? ? ? ? bfsize,p ? : ? Integer; ? //中間變量 ?
? ? ? ? ? ? ? ? ? angle ? : ? Double; ? //角度 ?
? begin ?
? ? ? ? ? ? ? ? ? count:= ? 1 ? shl ? r; ? //傅立葉變換點數 ?
? ? ? ? ? ? ? ? ? // ? 分配運算所需存儲器 ?
? ? ? ? ? ? ? ? ? SetLength(W, ? count ? div ? 2); ?
? ? ? ? ? ? ? ? ? SetLength(X1,count); ?
? ? ? ? ? ? ? ? ? SetLength(X2,count); ?
? ? ? ? ? ? ? ? ? // ? 計算加權系數 ?
? ? ? ? ? ? ? ? ? for ? i:=0 ? to ? count ? div ? 2 ? -1 ? do ?
? ? ? ? ? ? ? ? ? begin ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? angle ? := ? -i ? * ? PI ? * ? 2 ? / ? count; ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? W[i].vReal ? := ? cos(angle); ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? W[i].vImaginary ? := ? sin(angle); ?
? ? ? ? ? ? ? ? ? end; ?
? ? ? ? ? ? ? ? ? // ? 將時域點寫入X1 ?
? ? ? ? ? ? ? ? ? for ? i:=0 ? to ? count-1 ? do ?
? ? ? ? ? ? ? ? ? begin ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? X1[i].vReal:=TD[i]; ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? X1[i].vImaginary:=0; ?
? ? ? ? ? ? ? ? ? end; ?
? ? ? ? ? ? ? ? ? // ? 采用蝶形算法進行快速傅立葉變換 ? ?
? ? ? ? ? ? ? ? ? for ? k:=0 ? to ? r-1 ? do ?
? ? ? ? ? ? ? ? ? begin ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? for ? j:=0 ? to ? (1 ? shl ? k)-1 ? do ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? begin ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? bfsize:= ? 1 ? shl ? (r-k); ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? for ? i:=0 ? to ? bfsize ? div ? 2 ? -1 ? do ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? begin ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? p ? := ? j ? * ? bfsize; ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? X2[i+p]:=ComplexAdd(X1[i+p],X1[i+p+bfsize ? div ? 2]); ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? X2[i+p+bfsize ? div ? 2]:=ComplexMultiply( ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ComplexMinus(X1[i+p],X1[i+p+ ? bfsize ? div ? 2]), ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? W[i*(1 ? shl ? k)]); ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? end; ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? end; ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ComplexArrayWrap(X1,X2); ?
? ? ? ? ? ? ? ? ? end; ?
? ? ? ? ? ? ? ? ? // ? 重新排序 ? ?
? ? ? ? ? ? ? ? ? Setlength(Fd,count ? div ? 2); ?
? ? ? ? ? ? ? ? ? for ? j:=0 ? to ? count ? div ? 2 ? - ? 1 ? do ?
? ? ? ? ? ? ? ? ? begin ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? p ? := ? 0; ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? for ? i:=0 ? to ? r-1 ? do ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? begin ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? if ? (j ? and ? (1 ? shl ? i))>0 ? then ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? p:=p+(1 ? shl ? (r-i-1)); ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? end; ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? FD[j]:=X1[p].vReal; ?
? ? ? ? ? ? ? ? ? end; ?
? ? ? ? ? ? ? ? ? // ? 釋放內存 ? ?
? ? ? ? ? ? ? ? ? Setlength(W,0); ?
? ? ? ? ? ? ? ? ? Setlength(X1,0); ?
? ? ? ? ? ? ? ? ? Setlength(X2,0); ?
? end; ?
? 上面是fft的算法實現代碼,但是我在調用的時候,發現TD,FD輸入參數有些問題, ?
? 我原始數據是一個橋梁檢測的信號,一個分量是頻率,一個分量是幅值,我要如何通過FFT算法,把這個信號圖,分解成為頻率圖和幅度圖? ?
? 請教各位高手,謝謝!
http://www.delphi2007.net/DelphiBase/html/delphi_20061210124307209.html
type ? ? TComplexData ? = ? record ? //自定義復數類型 ?
? ? ? ? ? ? vReal,vImaginary: ? Double; ?
? end; ?
? TComplexDataArray ? = ? array ? of ? TComplexData; ? //實數數組 ? ?
? TExtendedArray ? = ? array ? of ? Extended; ? //復數數組 ?
? ?
? function ? ComplexAdd(c1,c2:TComplexData):TComplexData; ? ? ? ? ? ? ? //復數相加 ?
? function ? ComplexMinus(c1,c2:TComplexData):TComplexData; ? ? ? ? ? //復數相減 ?
? function ? ComplexMultiply(c1,c2:TComplexData):TComplexData; ? ? //復數相乘 ?
? procedure ? ComplexArrayWrap(var ? a1,a2:TComplexDataArray); ? ? ? ? //復數數組互換 ?
? procedure ? DoFFT(TD:TExtendedArray; ? var ? FD:TExtendedArray; ? r:Integer); ?
? ?
? procedure ? Tmainform.DoFFT(TD:TExtendedArray; ? var ? FD:TExtendedArray; ? r:Integer); ?
? var ?
? ? ? ? ? ? ? ? ? i,j,k ? : ? Integer; ? //循環變量 ?
? ? ? ? ? ? ? ? ? bfsize,p ? : ? Integer; ? //中間變量 ?
? ? ? ? ? ? ? ? ? angle ? : ? Double; ? //角度 ?
? begin ?
? ? ? ? ? ? ? ? ? count:= ? 1 ? shl ? r; ? //傅立葉變換點數 ?
? ? ? ? ? ? ? ? ? // ? 分配運算所需存儲器 ?
? ? ? ? ? ? ? ? ? SetLength(W, ? count ? div ? 2); ?
? ? ? ? ? ? ? ? ? SetLength(X1,count); ?
? ? ? ? ? ? ? ? ? SetLength(X2,count); ?
? ? ? ? ? ? ? ? ? // ? 計算加權系數 ?
? ? ? ? ? ? ? ? ? for ? i:=0 ? to ? count ? div ? 2 ? -1 ? do ?
? ? ? ? ? ? ? ? ? begin ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? angle ? := ? -i ? * ? PI ? * ? 2 ? / ? count; ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? W[i].vReal ? := ? cos(angle); ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? W[i].vImaginary ? := ? sin(angle); ?
? ? ? ? ? ? ? ? ? end; ?
? ? ? ? ? ? ? ? ? // ? 將時域點寫入X1 ?
? ? ? ? ? ? ? ? ? for ? i:=0 ? to ? count-1 ? do ?
? ? ? ? ? ? ? ? ? begin ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? X1[i].vReal:=TD[i]; ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? X1[i].vImaginary:=0; ?
? ? ? ? ? ? ? ? ? end; ?
? ? ? ? ? ? ? ? ? // ? 采用蝶形算法進行快速傅立葉變換 ? ?
? ? ? ? ? ? ? ? ? for ? k:=0 ? to ? r-1 ? do ?
? ? ? ? ? ? ? ? ? begin ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? for ? j:=0 ? to ? (1 ? shl ? k)-1 ? do ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? begin ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? bfsize:= ? 1 ? shl ? (r-k); ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? for ? i:=0 ? to ? bfsize ? div ? 2 ? -1 ? do ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? begin ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? p ? := ? j ? * ? bfsize; ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? X2[i+p]:=ComplexAdd(X1[i+p],X1[i+p+bfsize ? div ? 2]); ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? X2[i+p+bfsize ? div ? 2]:=ComplexMultiply( ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ComplexMinus(X1[i+p],X1[i+p+ ? bfsize ? div ? 2]), ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? W[i*(1 ? shl ? k)]); ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? end; ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? end; ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ComplexArrayWrap(X1,X2); ?
? ? ? ? ? ? ? ? ? end; ?
? ? ? ? ? ? ? ? ? // ? 重新排序 ? ?
? ? ? ? ? ? ? ? ? Setlength(Fd,count ? div ? 2); ?
? ? ? ? ? ? ? ? ? for ? j:=0 ? to ? count ? div ? 2 ? - ? 1 ? do ?
? ? ? ? ? ? ? ? ? begin ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? p ? := ? 0; ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? for ? i:=0 ? to ? r-1 ? do ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? begin ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? if ? (j ? and ? (1 ? shl ? i))>0 ? then ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? p:=p+(1 ? shl ? (r-i-1)); ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? end; ?
? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? FD[j]:=X1[p].vReal; ?
? ? ? ? ? ? ? ? ? end; ?
? ? ? ? ? ? ? ? ? // ? 釋放內存 ? ?
? ? ? ? ? ? ? ? ? Setlength(W,0); ?
? ? ? ? ? ? ? ? ? Setlength(X1,0); ?
? ? ? ? ? ? ? ? ? Setlength(X2,0); ?
? end; ?
? 上面是fft的算法實現代碼,但是我在調用的時候,發現TD,FD輸入參數有些問題, ?
? 我原始數據是一個橋梁檢測的信號,一個分量是頻率,一個分量是幅值,我要如何通過FFT算法,把這個信號圖,分解成為頻率圖和幅度圖? ?
? 請教各位高手,謝謝!
網上有很多快速傅立葉變換的例子. ? 有DELPHI實現的算法庫
有源碼庫的
噢 ? 來湊個熱鬧
現在的問題是不知道怎么初始化參數,也就是TD和FD,有沒有這方面的實例?謝謝! ?
? 我的Email是:chenamo5776@126.com
總結
- 上一篇: 产品经理在跨部门沟通中常见问题和解决办法
- 下一篇: 实施TDD时的常见问题