FFT插值算法程序.docx
《FFT插值算法程序.docx》由会员分享，可在线阅读，更多相关《FFT插值算法程序.docx（51页珍藏版）》请在冰豆网上搜索。
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.
升级会员 