FFT插值算法程序.docx
《FFT插值算法程序.docx》由会员分享,可在线阅读,更多相关《FFT插值算法程序.docx(51页珍藏版)》请在冰豆网上搜索。
![FFT插值算法程序.docx](https://file1.bdocx.com/fileroot1/2023-1/10/f0772132-5c4a-4aee-ad7f-a69b90768a96/f0772132-5c4a-4aee-ad7f-a69b90768a961.gif)
FFT插值算法程序
1.#include "stdafx.h"
2.#include "Fft.h"
3.#include "math.h"
4.
5.const double MOD_MAX = 65535.0;
6.const double DISP_MAX = 1.0/255.0;
7.extern FILE *fp;
8.
9.
10.
11.
12.
13.//FFT运算必须参数
14.int fft_point,fft_order,fft_divide,fft_window,fft_scale;
15.bool fft_cover;
16.float filter[7];//FIR滤波参数
17.
18.extern FILE *fpIandQ;
19.extern bool m_bIqWrite;
20.double prFilter[256],piFilter[256];
21.
22.
23./**************************************************************************************
24. 0 RectangleWindow矩形窗
25. FFT变换结果为对称型,矩形窗是使信号突然截断,旁瓣会很大,且衰减较慢,旁瓣的第一个负
26. 峰值为主瓣的21%,第一个正峰值为主瓣的12.6%,第二个负峰值为主瓣的9%,效果一般,泄
27. 漏较大。
28.**************************************************************************************/
29.double WINAPI RectangleWindow(int t)
30.{
31. double wt;
32. wt=1.0;
33. return wt;
34.}
35.
36./**************************************************************************************
37. 1 TriangleWindow三角窗,也称费杰(Fejer)窗,Bartlett
38.**************************************************************************************/
39.double WINAPI TriangleWindow(int t)
40.{
41.
42. double wt;
43. wt=1-t/fft_point;
44. return wt;
45.}
46.
47./**************************************************************************************
48. 2 HanningWindwo汉宁窗,即升余弦窗
49.**************************************************************************************/
50.double WINAPI HanningWindow(int t)
51.{
52. double wt;
53. wt=(1-cos(2*PI*t/fft_point))/2;
54. return wt;
55.}
56.
57./**************************************************************************************
58. 3 HammingWindow海明窗,即改进的升余弦窗
59.**************************************************************************************/
60.double WINAPI HammingWindow(int t)
61.{
62. double wt;
63. wt=0.54-0.46*cos(2*PI*t/fft_point);
64. return wt;
65.}
66.
67./**************************************************************************************
68. 4 BlackmanWindow布来克曼窗,即三阶升余弦窗
69.**************************************************************************************/
70.double WINAPI BlackmanWindow(int t)
71.{
72. double wt;
73. wt=0.42-0.5*cos(2*PI*t/fft_point)+0.08*cos(4*PI*t/fft_point);
74. return wt;
75.}
76.
77./**************************************************************************************
78. 5 CosgradeWindow余弦坡度窗
79.**************************************************************************************/
80.double WINAPI CosgradeWindow(int t)
81.{
82. double wt;
83. if(t= int(4*fft_point/5) )
84. wt=1.0;
85. else
86. wt=(1+cos(5*PI*t/fft_point))/2;
87. return wt;
88.}
89.
90.
91./**************************************************************************************
92. 6 ParzenWindow帕曾窗
93.**************************************************************************************/
94.double WINAPI ParzenWindow(int t)
95.{
96. double wt;
97. if(t= int(fft_point/2) )
98. wt=1-6*pow(t/fft_point,2)+6*pow(t/fft_point,3);
99. else
100. wt=2*pow(1-t/fft_point,3);
101. return wt;
102.}
103.
104./**************************************************************************************
105. 7 ExponentWindow指数窗
106.**************************************************************************************/
107.double WINAPI ExponentWindow(int t)
108.{
109. double wt,alf;
110. alf=0.0001;
111. wt=exp((-1.0)*alf*t);
112. return wt;
113.}
114.
115./**************************************************************************************
116. 8 GaussWindow高斯窗
117.**************************************************************************************/
118.double WINAPI GaussWindow(int t)
119.{
120. double wt,alf;
121. alf=0.0001;
122. wt=exp(-1*alf*t*t);
123. return wt;
124.}
125.
126./**************************************************************************************
127. 9 KaiserWindow凯泽窗,
128.凯泽在1966(1974)发现,利用第一类零阶修正(变形)贝赛尔函数可以构成一种近似最佳的窗
129.最主要参数:
beat-可同时调整主瓣宽度和旁瓣,beat越大,窗越窄,频谱旁瓣越小,但主瓣相应增加
130. beta=0相当与矩形窗
131.**************************************************************************************/
132.double WINAPI KaiserWindow(int t)
133.{
134. double wt,alfa,beta;
135. alfa=(fft_point-1)/2;
136. beta=0;
137. //Bessel零阶第一类贝塞尔函数
138. wt=Jim_Bessel0_R(0,beta*sqrt(1-pow((t-alfa)/alfa,2)))/Jim_Bessel0_R(0,beta);
139. return wt;
140.}
141.
142./**************************************************************************************
143. 10 ChebWindow契比雪夫窗
144.**************************************************************************************/
145.double WINAPI ChebWindow(int t)
146.{
147. double wt;
148. wt=1.0;
149. return wt;
150.}
151.
152./**************************************************************************************
153. 11 BartlettWindow巴特利特窗
154.**************************************************************************************/
155.double WINAPI BartlettWindow(int t)
156.{
157. double wt;
158. wt=1.0;
159. return wt;
160.}
161.
162.bool WINAPI fftInit(int point,int order,int divide,bool cover,int scale)
163.{
164.
165. fft_point=point;
166. fft_order=order;
167. fft_divide=divide;
168. fft_cover=cover;
169. fft_scale=scale;
170. //采用8阶的FIR滤波
171. Jim_FirFilter(7,3,50,7800/2.0,7800.0,1,filter);
172. return TRUE;
173.}
174.
175.
176.BOOL WINAPI fftEnd()
177.{
178. return TRUE;
179.}
180.
181.
182.
183.
184.//傅立叶变换
185.BOOL WINAPI fftTransform(VOID *inBuf,VOID *outBuf,int windows,int scale,int nProbe,int nWork)
186.{
187. BYTE *srcPtr = (BYTE*)inBuf;//数据源1024字节
188. BYTE *dstPtr = (BYTE*)outBuf;//数据果256字节
189. //double *testPtr =(double *)testBuf;
190. WORD Ivalue,Qvalue;
191. //WORD IvalueRev,QvalueRev;
192. unsigned char flagQ,flagI;
193. //double alfa,e;
194. int i,j,Iorg=0,Qorg=0;
195. double mod = 0;
196. if(nProbe==LEFT_PROBE)
197. {
198. flagQ=0x30,flagI=0x20;
199. }
200. else
201. {
202. flagQ=0x10,flagI=0x0;
203. }
204.
205. //每次取4个字节数据,分离I/Q分量,判断db(12),0表示I分量,1表示Q分量
206. if (((srcPtr[1]&0xf0) ==flagQ) && ((srcPtr[3]&0xf0) ==flagI))
207. {
208. Iorg = 2;
209. Qorg = 0;
210. }
211. if (((srcPtr[1]&0xf0) ==flagI) && ((srcPtr[3]&0xf0) ==flagQ))
212. {
213. Iorg = 0;
214. Qorg = 2;
215. }
216. double pr[256],pi[256],w;
217.
218.
219. //计算自相关函数和互相关函数
220. //double Rii=0.0,Rqq=0.0,Riq=0.0;
221.
222. for ( i=0;i223. {
224.
225. //更改存储的偏移地址分离I/Q,
226. Ivalue = *((short*)(srcPtr + i*4 + Iorg));
227. Qvalue = *((short*)(srcPtr + i*4 + Qorg));
228. //~位非运算,低12位位有效数据
229. pr[i]=double(Ivalue & 0xfff);
230. pi[i]=double(Qvalue & 0xfff);
231.
232. /*Rii=Rii+double(pr[i]*pr[i]);
233. Rqq=Rqq+double(pi[i]*pi[i]);
234. Riq=Riq+double(pr[i]*pi[i]);*/
235. }
236.
237. /*Rii=Rii/fft_point;
238. Rqq=Rqq/fft_point;
239. Riq=Riq/fft_point;
240.
241. alfa=asin(Riq/sqrt(Rii*Rqq));
242. e=sqrt(Rqq/Rii)-1;*/
243.
244. for ( i=0;i245. {
246. //防止谱泄漏,进行加窗处理
247. switch(windows)
248. {
249. case WND_RECTANGLE:
w=RectangleWindow(i); break;
250. case WND_TRIANGLE:
w=TriangleWindow(i); break;
251. case WND_HANNING:
w=HanningWindow(i); break;
252. case WND_HAMMING:
w=HammingWindow(i); break;
253. case WND_BLACKMAN:
w=BlackmanWindow(i); break;
254. case WND_COSGRADE:
w=CosgradeWindow(i); break;
255. case WND_PARZEN:
w=ParzenWindow(i); break;
256. case WND_EXPONENT:
w=ExponentWindow(i); break;
257. case WND_GAUSS:
w=GaussWindow(i); break;
258. case WND_KAISER:
w=KaiserWindow(i); break;
259. case WND_CHEB:
w=ChebWindow(i); break;
260. case WND_BARTLETT:
w=BartlettWindow(i); break;
261. default:
w=RectangleWindow(i); break;
262. }
263.
264. /*pr[i]=((1+e)*cos(alfa)*pr[i])/((1+e)*cos(alfa))*w;
265. pi[i]=(pi[i]-(1+e)*sin(alfa)*pr[i])/((1+e)*cos(alfa))*w;*/
266.
267. pr[i]=pr[i]*w;
268. pi[i]=pi[i]*w;
269. prFilter[i]=piFilter[i]=0.0;
270. for(j=0;j<7;j++)
271. {
272. if(i+j>255)
273. {
274.