开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 数学计算 > 查看源码
SignalProcess.cpp
资源名称：SignalProcess.rar　[点击查看]
上传用户：hldfqc
上传日期：2022-07-25
资源大小：12k
文件大小：10k
源码类别：

数学计算

开发平台：

Visual C++

SignalProcess.cpp：源码内容
  1. // SignalProcess.cpp : Defines the entry point for the DLL application.
  2. //
  4. #include "stdafx.h"
  5. #include "SignalProcess.h"
  7. BOOL APIENTRY DllMain( HANDLE hModule, 
  8.                        DWORD  ul_reason_for_call, 
  9.                        LPVOID lpReserved
  10.  )
  11. {
  12.     return TRUE;
  13. }
  14. /*-----------------------------------------------  
  15.            FFT函数
  16.     data:指向数据序列地指针
  17. a   :指向data的DFT序列的指针
  18. L   :2的L次方为FFT的点数
  19. --------------------------------------------------*/
  20. int fft(int *data,complex <double> *a,int L)
  21. {
  22. complex <double> u;
  23. complex <double> w;
  24. complex <double> t;
  25. unsigned n=1,nv2,nm1,k,le,lei,ip;
  26. int i,j,m,number;
  27. double tmp;
  28. n<<=L;
  29. for(number = 0; number<n; number++)
  30. {
  31. a[number] = complex <double> (data[number],0);
  32. }
  33. nv2=n>>1;
  34. nm1=n-1;
  35. j=0;
  36. for(i=0;i<nm1;i++)
  37. {
  38. if(i<j)
  39. {
  40. t=a[j];
  41. a[j]=a[i];
  42. a[i]=t;
  43. }
  44. k=nv2;
  45. while(k<=j)
  46. {
  47. j-=k;
  48. k>>=1;
  49. }
  50. j+=k;
  51. }
  52. le=1;
  53. for(m=1;m<=L;m++)
  54. {
  55. lei=le;
  56. le<<=1;
  57. u=complex<double>(1,0);
  58. tmp=PI/lei;
  59. w=complex<double>(cos(tmp),-sin(tmp));
  60. for(j=0;j<lei;j++)
  61. {
  62. for(i=j;i<n;i+=le)
  63. {
  64. ip=i+lei;
  65. t=a[ip]*u;
  66. a[ip]=a[i]-t;
  67. a[i]+=t;
  68. }
  69. u*=w;
  70. }
  71. }
  73. return 0;
  74. }
  75. /*-----------------------------------------------  
  76.            IFFT函数
  77.     data:指向数据序列地指针
  78. a   :指向data的DFT序列的指针
  79. L   :2的L次方为FFT的点数
  80. --------------------------------------------------*/
  82. int ifft(complex <double> *a,int *data,int L)
  83. {
  84. complex <double> u;
  85. complex <double> w;
  86. complex <double> t;
  87. unsigned n=1,nv2,nm1,k,le,lei,ip;
  88. int i,j,m,number;
  89. double tmp;
  90. n<<=L;
  91. nv2=n>>1;
  92. nm1=n-1;
  93. j=0;
  94. for(i=0;i<nm1;i++)
  95. {
  96. if(i<j)
  97. {
  98. t=a[j];
  99. a[j]=a[i];
  100. a[i]=t;
  101. }
  102. k=nv2;
  103. while(k<=j)
  104. {
  105. j-=k;
  106. k>>=1;
  107. }
  108. j+=k;
  109. }
  110. le=1;
  111. for(m=1;m<=L;m++)
  112. {
  113. lei=le;
  114. le<<=1;
  115. u=complex<double>(1,0);
  116. tmp=PI/lei;
  117. w=complex<double>(cos(tmp),sin(tmp));
  118. for(j=0;j<lei;j++)
  119. {
  120. for(i=j;i<n;i+=le)
  121. {
  122. ip=i+lei;
  123. t=a[ip]*u;
  124. a[ip]=a[i]-t;
  125. a[i]+=t;
  126. }
  127. u*=w;
  128. }
  129. }
  130. for(number = 0; number<n; number++)
  131. {
  132. data[number] = ceil(a[number].real())/n;
  133. a[number] = a[number]/complex<double>(n,0);
  134. }
  135. return 0;
  136. }
  138. /*--------------------------------------------------------------
  139.                      Hilbert变换函数
  140. data：指向信号序列的指针
  141. filterdata：指向包络序列的指针
  142. dn：信号序列的点数
  143. ----------------------------------------------------------------*/
  144. int  __stdcall hilbert(int * data , int *filterdata,int dn)
  145. {
  146. int i = 0, j = 0, k = 0,l = 0,N = 0;
  147. complex<double> *zk;
  148. int *ldata;
  149. l = (int)(log(dn)/log(2))+1;
  150. N =(int) pow(2,l);
  151. zk = (complex<double>*)malloc(N*sizeof(complex<double>));
  152. ldata = (int *)malloc(N*sizeof(int));
  153. memcpy(ldata,data,dn*sizeof(int));
  154. for(i=dn;i<N;i++)
  155. {
  156. ldata[i] = 0;
  157. }
  158. fft(ldata,zk,l);
  159. for(i=0;i<N;i++)
  160. {
  161. if(i>=1 && i<=N/2-1)
  162. {
  163. zk[i] = complex<double>(2,0)*zk[i];
  164. }
  165. if(i>=N/2 && i<=N-1)
  166. {
  167. zk[i]= complex<double> (0,0);
  168. }
  169. }
  170. ifft(zk,ldata,l);
  171.     for(i = 0 ;i<dn;i++)
  172. {
  173. filterdata[i] = (int)sqrt(pow(zk[i].imag(),2)+pow(zk[i].real(),2));
  174. }
  175. free(zk);
  176. free(ldata);
  177. return 0;
  178. }
  180. int __stdcall conv(int *h,int *data,int *result,int hn,int dn)
  181. {
  182. int l,i,j,k,N;
  183. complex<double> *hk,*datak;
  184. l = (int)(log(hn+dn-1)/log(2))+1;
  185. N =(int) pow(2,l);
  186. int *lh,*ldata;
  187. lh =(int*)malloc(N*sizeof(int));
  188. ldata =(int*)malloc(N*sizeof(int));
  189. memcpy(lh,h,hn*sizeof(int));
  190. memcpy(ldata,data,dn*sizeof(int));
  191. for(i=hn;i<N;i++)
  192. {
  193. lh[i] = 0;
  194. }
  195. for(j=dn;j<N;j++)
  196. {
  197. ldata[j] = 0;
  198. }
  199.     hk = (complex <double> *) malloc(N*sizeof(complex<double>));
  200.     datak = (complex <double> *) malloc(N*sizeof(complex<double>));
  201. fft(lh,hk,l);
  202. fft(ldata,datak,l);
  203. for(k=0;k<N;k++)
  204. {
  205. datak[k] = datak[k]*hk[k];
  206. }
  207. ifft(datak,result,l);
  209. free(lh);
  210. free(ldata);
  211. free(hk);
  212. free(datak);
  213. return 0;
  214. }
  215. //////////////////////////////////////////////////////////////////
  216. int fft_f(double *data,complex <double> *a,int L)
  217. {
  218. complex <double> u;
  219. complex <double> w;
  220. complex <double> t;
  221. unsigned n=1,nv2,nm1,k,le,lei,ip;
  222. int i,j,m,number;
  223. double tmp;
  224. n<<=L;
  225. for(number = 0; number<n; number++)
  226. {
  227. a[number] = complex <double> (data[number],0);
  228. }
  229. nv2=n>>1;
  230. nm1=n-1;
  231. j=0;
  232. for(i=0;i<nm1;i++)
  233. {
  234. if(i<j)
  235. {
  236. t=a[j];
  237. a[j]=a[i];
  238. a[i]=t;
  239. }
  240. k=nv2;
  241. while(k<=j)
  242. {
  243. j-=k;
  244. k>>=1;
  245. }
  246. j+=k;
  247. }
  248. le=1;
  249. for(m=1;m<=L;m++)
  250. {
  251. lei=le;
  252. le<<=1;
  253. u=complex<double>(1,0);
  254. tmp=PI/lei;
  255. w=complex<double>(cos(tmp),-sin(tmp));
  256. for(j=0;j<lei;j++)
  257. {
  258. for(i=j;i<n;i+=le)
  259. {
  260. ip=i+lei;
  261. t=a[ip]*u;
  262. a[ip]=a[i]-t;
  263. a[i]+=t;
  264. }
  265. u*=w;
  266. }
  267. }
  269. return 0;
  270. }
  273. int conv_f(double *h,int *data,int *result,int hn,int dn)
  274. {
  275. int l,i,j,k,N;
  276. complex<double> *hk,*datak;
  277. l = (int)(log(hn+dn-1)/log(2))+1;
  278. N =(int) pow(2,l);
  279. int *ldata;
  280. double *lh;
  281. lh =(double*)malloc(N*sizeof(double));
  282. ldata =(int*)malloc(N*sizeof(int));
  283. memcpy(lh,h,hn*sizeof(double));
  284. memcpy(ldata,data,dn*sizeof(int));
  285. for(i=hn;i<N;i++)
  286. {
  287. lh[i] = 0;
  288. }
  289. for(j=dn;j<N;j++)
  290. {
  291. ldata[j] = 0;
  292. }
  293.     hk = (complex <double> *) malloc(N*sizeof(complex<double>));
  294.     datak = (complex <double> *) malloc(N*sizeof(complex<double>));
  295. fft_f(lh,hk,l);
  296. fft(ldata,datak,l);
  297. for(k=0;k<N;k++)
  298. {
  299. datak[k] = datak[k]*hk[k];
  300. }
  301. ifft(datak,result,l);
  302. free(lh);
  303. free(ldata);
  304. free(hk);
  305. free(datak);
  306. return 0;
  307. }
  309. /*int __stdcall firwin_e(int n,int band,int fl,int fh,int fs,int wn,int *h)
  310. {
  311. int i,n2,mid;
  312. double s,wc1,wc2,beta,delay;
  313. beta=0.0;
  314. double fln = (double)fl / fs;
  315. double fhn = (double)fh / fs;
  316. // if(wn==7)
  317. // {
  318. // printf("input beta parameter of Kaiser window(3<beta<10)n");
  319. // scanf("%lf",&beta);
  320. // }
  321. beta = 6;
  322. if((n%2)==0)
  323. {
  324. n2=n/2-1;
  325. mid=1;
  326. }
  327. else
  328. {
  329. n2=n/2;
  330. mid=0;
  331. }
  332. delay=n/2.0;
  333. wc1=2.0*PI*fln;
  334. if(band>=3) wc2=2.0*PI*fhn;
  335. switch(band)
  336. {
  337. case 1://低通
  338. {
  339. for (i=0;i<=n2;i++)
  340. {
  341. s=i-delay;
  342. *(h+i)=(int)((sin(wc1*s)/(PI*s))*window(wn,n+1,i,beta));
  343. *(h+n-i)=*(h+i);
  344. }
  345. if(mid==1) *(h+n/2)=(int)(wc1/PI);
  346. break;
  347. }
  349. case 2: //高通
  350. {
  351. for (i=0;i<=n2;i++)
  352. {
  353. s=i-delay;
  354. *(h+i)=(int)((sin(PI*s)-sin(wc1*s))/(PI*s));
  355. *(h+i)=*(h+i)*(int)(window(wn,n+1,i,beta));
  356. *(h+n-i)=*(h+i);
  357. }
  358. if(mid==1) *(h+n/2)=(int)(1.0-wc1/PI);
  359. break;
  360. }
  361. case 3: //带通
  362. {
  363. for (i=0;i<=n2;i++)
  364. {
  365. s=i-delay;
  366. *(h+i)=(int)((sin(wc2*s)-sin(wc1*s))/(PI*s));
  367. *(h+i)=*(h+i)*(int)(window(wn,n+1,i,beta));
  368. *(h+n-i)=*(h+i);
  369. }
  370. if(mid==1) *(h+n/2)=(int)((wc2-wc1)/PI);
  371. break;
  372. }
  373. case 4: //带阻
  374. {
  375. for (i=0;i<=n2;i++)
  376. {
  377. s=i-delay;
  378. *(h+i)=(int)((sin(wc1*s)+sin(PI*s)-sin(wc2*s))/(PI*s));
  379. *(h+i)=*(h+i)*(int)(window(wn,n+1,i,beta));
  380. *(h+n-i)=*(h+i);
  381. }
  382. if(mid==1) *(h+n/2)=(int)((wc1+PI-wc2)/PI);
  383. break;
  384. }
  385. }
  386. return 0;
  387. }*/
  389. double window(int type,int n,int i,double beta)
  390. {
  391. int k;
  392. double w=1.0;
  393. switch(type)
  394. {
  395. case 1:
  396. {
  397. w=1.0;
  398. break;
  399. }
  400. case 2:
  401. {
  402. k=(n-2)/10;
  403. if(i<=k) w=0.5*(1.0-cos(i*PI/(k+1)));
  404. if(i>n-k-2) w=0.5*(1.0-cos((n-i-1)*PI/(k+1)));
  405. break;
  406. }
  407. case 3:
  408. {
  409. w=1.0-fabs(1.0-2*i/(n-1.0));
  410. break;
  411. }
  412. case 4:
  413. {
  414. w=0.5*(1.0-cos(2*i*PI/(n-1)));
  415. break;
  416. }
  417. case 5:
  418. {
  419. w=0.54-0.46*cos(2*i*PI/(n-1));
  420. break;
  421. }
  422. case 6:
  423. {
  424. w=0.42-0.5*cos(2*i*PI/(n-1))+0.08*cos(4*i*PI/(n-1));
  425. break;
  426. }
  427. case 7:
  428. {
  429. w=kaiser(i,n,beta);
  430. break;
  431. }
  432. }
  433. return(w);
  434. }
  436. double kaiser(int i,int n,double beta)
  437. {
  438. double a,w,a2,b1,b2,beta1;
  439. b1=bessel0(beta);
  440. a=2.0*i/(double)(n-1)-1.0;
  441. a2=a*a;
  442. beta1=beta*sqrt(1.0-a2);
  443. b2=bessel0(beta1);
  444. w=b2/b1;
  445. return(w);
  446. }
  448. double bessel0(double x)
  449. {
  450. int i;
  451. double d,y,d2,sum;
  452. y=x/2.0;
  453. d=1.0;
  454. for(i=1;i<=25;i++)
  455. {
  456. d=d*y/i;
  457. d2=d*d;
  458. sum=sum+d2;
  459. if(d2<sum*(1.0e-8)) break;
  460. }
  461. return(sum);
  462. }
  465. int __stdcall firwin_e(int n,int band,int fl,int fh,int fs,int wn, int *data,int *result,int dn)
  466. {
  467. int i,n2,mid;
  468. double *h;
  469. double s,wc1,wc2,beta,delay;
  470. beta=0.0;
  471. double fln = (double)fl / fs;
  472. double fhn = (double)fh / fs;
  473. // if(wn==7)
  474. // {
  475. // printf("input beta parameter of Kaiser window(3<beta<10)n");
  476. // scanf("%lf",&beta);
  477. // }
  478. h = (double *)malloc((n+1)*sizeof(double));
  479. beta = 6;
  480. if((n%2)==0)
  481. {
  482. n2=n/2-1;
  483. mid=1;
  484. }
  485. else
  486. {
  487. n2=n/2;
  488. mid=0;
  489. }
  490. delay=n/2.0;
  491. wc1=2.0*PI*fln;
  492. // FILE *fp;
  493. if(band>=3) wc2=2.0*PI*fhn;
  494. switch(band)
  495. {
  496. case 1://低通
  497. {
  498. for (i=0;i<=n2;i++)
  499. {
  500. s=i-delay;
  501. *(h+i)=(sin(wc1*s)/(PI*s))*window(wn,n+1,i,beta);
  502. *(h+n-i)=*(h+i);
  503. }
  504. if(mid==1) *(h+n/2)=wc1/PI;
  505. break;
  506. }
  508. case 2: //高通
  509. {
  510. for (i=0;i<=n2;i++)
  511. {
  512. s=i-delay;
  513. *(h+i)=(sin(PI*s)-sin(wc1*s))/(PI*s);
  514. *(h+i)=*(h+i)*window(wn,n+1,i,beta);
  515. *(h+n-i)=*(h+i);
  516. }
  517. if(mid==1) *(h+n/2)=1.0-wc1/PI;
  518. break;
  519. }
  520. case 3: //带通
  521. {
  522. for (i=0;i<=n2;i++)
  523. {
  524. s=i-delay;
  525. *(h+i)=(sin(wc2*s)-sin(wc1*s))/(PI*s);
  526. *(h+i)=*(h+i)*window(wn,n+1,i,beta);
  527. *(h+n-i)=*(h+i);
  528. }
  529. if(mid==1) *(h+n/2)=(wc2-wc1)/PI;
  530. break;
  531. }
  532. case 4: //带阻
  533. {
  534. for (i=0;i<=n2;i++)
  535. {
  536. s=i-delay;
  537. *(h+i)=(sin(wc1*s)+sin(PI*s)-sin(wc2*s))/(PI*s);
  538. *(h+i)=*(h+i)*window(wn,n+1,i,beta);
  539. *(h+n-i)=*(h+i);
  540. }
  541. if(mid==1) *(h+n/2)=(wc1+PI-wc2)/PI;
  542. break;
  543. }
  544. }
  545. // fp=fopen("h.dat","w");
  546. // for(int p= 0;p<n+1;p++)
  547. // {
  548. // fprintf(fp, "%fn", h[p]);
  549. // }
  550. // fclose(fp);
  551. conv_f(h,data,result,n+1,dn);
  552. free(h);
  553. return 0;
  554. }