fft傅里叶变换问题
时间：10-02
整理：3721RD
点击：
现在固定频率8k采样，如果是50hz的信号相对应就是160个采样点，现在要进行傅里叶变换，有没有大神能帮我看看代码 pwcz函数是在干嘛，整个代码能不能帮我屡屡思路
- #include <stdio.h>
- #include <stdlib.h>
- #include <math.h>
-
- #define SAMP 160
- #define frequency 47.5 //为方便仿真频率先固定
- #define MN (1.0*SAMP*50/47.5)
- float tmp[21];
- float ratio[20];
- struct compx
- {
- float real;
- float imag;
- };
- struct compx sum[SAMP+1];
- float base; /*-基波-*/
- float component_O[20]; /* -各次谐波分量-*/
- float component_O_rt[20];
- float ftmp;
- float angle[20];
- int N = SAMP;
- #define RN SAMP*2
- float AD_RESULT[RN];
- float CZ_RESULT[RN];
-
- /*float amend_mt[20]= //谐波修正系数
- {
- 1.0,0.9752,0.9752,0.9752,0.9752,0.9752,0.9752,0.9752,0.9752,0.9752,
- 0.9752*1.08,0.9752*1.015,0.9752*1.022,0.9752*1.025,0.9752*1.030,
- 0.9752*1.032,0.9752*1.042,0.9752*1.055,0.9752*1.065,0.9752*1.083
- };
- */
-
-
-
- struct compx EE(struct compx b1,struct compx b2)
- {
- struct compx b3;
- b3.real=(b1.real*b2.real-b1.imag*b2.imag);
- b3.imag=(b1.real*b2.imag+b1.imag*b2.real);
- return(b3);
- }
-
- float pwcz(float x)
- {
- float y;
- float e1,e2,e3;
- unsigned int k;
-
- if((unsigned int)(x*10)%10<=5)
- {
- k=(unsigned int)x;
- printf("k=%d\t",k);
- }
- else
- {
- k=(unsigned int)x+1;
- printf("k1=%d\t",k);
- }
- e1=(x-k)*(x-(k+1))/2.0;
- e2=(x-(k-1))*(x-(k+1))/-1.0;
- e3=(x-(k-1))*(x-k)/2.0;
- y=e1*AD_RESULT[k-1]+e2*AD_RESULT[k]+e3*AD_RESULT[k+1];
- printf("y=%f\n", y);
- return y;
- }
-
- static void draw_data(float *buf,int len,int mode)
- {
- FILE *stream;
- float u;
- float temp;
- stream = fopen( "data.txt", "w" );
- int n = 0;
- while(len> 0)
- {
- temp = *buf;
- if( temp < 0)
- {
- temp = temp/4 + 60;
- }
- else
- {
- temp = temp/4 + 60;
- }
-
- fprintf(stream, "%d %f\t", n++, *buf);
- for (u = 0; u < temp-1; u++)
- {
- fprintf( stream, " ");
- }
- fprintf( stream, "*\n");
- buf += 1;
- len -= 1;
- };
- fclose( stream );
- }
-
-
-
- void main()
- {
- for(int n=1;n<=2*N;n++) // 正弦波加谐波计算测试
- {
- AD_RESULT[n]=220*sin(2*3.1415926*(n-1)/MN)
- +0*sin(2*2*3.1415926*(n-1)/MN) //2次
- // +22*sin(3*2*3.1415926*(n-1)/MN) //3次
- // +22*sin(4*2*3.1415926*(n-1)/MN) //4次
- // +22*sin(5*2*3.1415926*(n-1)/MN) //5次
- // +22*sin(6*2*3.1415926*(n-1)/MN) //6次
- // +22*sin(7*2*3.1415926*(n-1)/MN) //5次
- // +22*sin(8*2*3.1415926*(n-1)/MN) //5次
- // +22*sin(9*2*3.1415926*(n-1)/MN) //5次
- // +6.6*sin(10*2*3.1415926*(n-1)/MN) //10次
- // +6.6*sin(11*2*3.1415926*(n-1)/MN) //5次
- // +6.6*sin(12*2*3.1415926*(n-1)/MN) //5次
- // +6.6*sin(13*2*3.1415926*(n-1)/MN) //5次
- // +6.6*sin(14*2*3.1415926*(n-1)/MN) //5次
- // +6.6*sin(15*2*3.1415926*(n-1)/MN) //5次
- // +6.6*sin(16*2*3.1415926*(n-1)/MN) //16次
- +11*sin(17*2*3.1415926*(n-1)/MN)
- +11*sin(18*2*3.1415926*(n-1)/MN)
- +11*sin(21*2*3.1415926*(n-1)/MN);
- // +22*sin(19*2*3.1415926*(n-1)/MN); //19次
- // printf("AD_RESULT[%d]=%f\n", n, AD_RESULT[n]);
-
-
- }
- for(n=1;n<=SAMP;n++) // 正弦波加谐波计算测试
- {
- printf("AD_RESULT[%d]=%f\t", n, AD_RESULT[n]);
- CZ_RESULT[n]=pwcz(n*50.0/frequency);
- }
-
- float tmp2;
- float temp;
-
- for(int m=1; m<22;m++)
- {
- for(n=1; n<=SAMP; n++)
- {
- sum[m].real+=CZ_RESULT[n]*cos(2*3.1415926*m*n/SAMP);
- }
- for(n=1; n<=SAMP; n++)
- {
- sum[m].imag+=CZ_RESULT[n]*sin(2*3.1415926*m*n/SAMP);
- }
-
- tmp2=2*sqrt((sum[m].real*sum[m].real+sum[m].imag*sum[m].imag))/SAMP;
- temp = tmp2;
- component_O[m]=temp;
- ftmp = component_O[m] * component_O[m];
- component_O_rt[m] = sqrt(ftmp);
- ratio[m] = component_O_rt[m]*100/component_O_rt[1];
- // printf("component_O_rt[%d]=%f ratio[%d] = %f angle[%d]=%f\n", n, component_O_rt[n], n, ratio[n], n, angle[n]);
- printf("component_O_rt[%d]=%f ratio[%d] = %f \n", m, component_O_rt[m], m, ratio[m]);
- }
- draw_data(AD_RESULT, 2*SAMP, 0);
- draw_data(CZ_RESULT, SAMP, 0);
- }