parametet.m
clc; clear; load('src.mat') CZT_N = 2048; CZT_M = 2048; CZT_W = exp(-j*(2*pi/8192)); CZT_A = exp(j*1*pi/256);
该文件罗列了用到的参数的取值情况，src.mat文件是matlab保存下来的原始数据，数据情况：2048*1个数据，每个数据都是科学计数法表示的。
load命令，可以将文件中的数据加载到项目中，以便后续使用（如何使用，还依赖一个模块化仿真图）。
czt_w/A都是指数形式的复数：在c++语言中复数运算库——<complex>中定义了复数的运算，其中包括复数的指数运算——pow(a,b)=a^b和exp(a)=e^a。
c++中定义复数的形式：complex <double> var_name(real_port,imag_port);
这类似于在定义一个复数变量时调用了一个给实部和虚部赋值的函数，其实在c++中复数是一个类的定义，定义复数就是定义对象，给复数赋值也就只能是初始化方式赋值了——调用构造函数赋值，第二种赋值方式就是对象赋值。（目前了解的。。。。）
定义复数类的形式（形似）
class complex
{
public:
double real_port;
double imag_port;
};
例如：
complex <double> abc(1.2,3.12);
那么：
abc=(1.2,3.12)
real(abc)=1.2; imag(abc)=3.12
complex <double> bbb; bbb=abc; -->bbb=(1.2,3.12)
有了这点理解，所以：
parameter.h就好写了
1 #include <complex> 2 #include <vector> 3 #include <iostream> 4 5 using namespace std; 6 7 double pi = 3.14159265359; 8 9 int czt_n = 2048; 10 int czt_m = 2048; 11 complex <double> czt_w_1(0,-2*pi/8192); 12 complex <double> czt_w = exp(czt_w_1); 13 complex <double> czt_a_1(0,pi/256); 14 complex <double> czt_a = exp(czt_a_1);
CZT.m:
function y=CZT(u,CZT_N,CZT_M,CZT_W,CZT_A) L = CZT_M + CZT_N ; %L=4096 g = zeros(L,1); %g为L行1列的全0矩阵 for i=0:CZT_N-1 %矩阵g元素赋值 g(i+1,1) = CZT_A^(-i)*CZT_W^(i^2/2)*u(i+1,1); end G=fft(g); h = zeros(L,1); for i=0:CZT_M-1 h(i+1,1)=CZT_W^(-i^2/2); end for i=L-CZT_N+1:L-1 h(i+1,1)=CZT_W^(-(L-i)^2/2); end H=fft(h); Q=H.*G; q=ifft(Q); X=zeros(CZT_M,1); for i=0:CZT_M-1 X(i+1,1)=CZT_W^(i^2/2)*q(i+1); end y=X;
czt.m文件定义了一个函数——czt。（文件名与函数名相同）
function是matlab中定义函数的关键字，y表示返回值，u,CZT_N,CZT_M,CZT_W,CZT_A，表示参数，那个这个函数的理解就是用：u,CZT_N,CZT_M,CZT_W,CZT_A这些参数经过一系列步骤得到输出y的过程，这些步骤就是函数体，其中u表示什么呢？，它就是一个矩阵，元素为src.mat文件中的数据，因此u是一个2048*1的矩阵（矩阵可以用数组在c语言环境下代替）。
统领全篇：用到的计算方法有：fft、ifft、矩阵的点乘、复数的指数运算。
L = CZT_M + CZT_N ; %L=4096 g = zeros(L,1); %g为L行1列的全0矩阵 for i=0:CZT_N-1 %矩阵g元素赋值 g(i+1,1) = CZT_A^(-i)*CZT_W^(i^2/2)*u(i+1,1); end G=fft(g);
在写这部分代码时，首先要做的几个准备工作：
1、g矩阵是一个4096*1的矩阵，可以用一维数组代替，数组的元素性质是复数。
2、u矩阵的元素在这之前还是以文件的形式存储的，要读到u变量中去（src_real.txt只记录了实部，src_imag.txt只记录了虚部）
3、g = zeros(L,1)表示的是变量初始化过程，因此数组还得初始化。
4、值得注意的是这里的G由于要在下面参与运算，所以不能是函数的局部变量，还能提取fft函数，使其能在后续计算中能被访问。
1 // test_9.cpp : 定义控制台应用程序的入口点。 2 // 3 4 #include "stdafx.h" 5 #include <vector> 6 #include <complex> 7 #include <iostream> 8 #include "parameter.h" 9 #include "fftw3.h" 10 #include <stdlib.h> 11 #include <malloc.h> 12 13 14 complex <double> u[2048]; 15 16 /***********读原始数据到矩阵（数组u）***************/ 17 void src_txt_u() 18 { 19 int i=0; 20 double u_real=0; 21 double u_imag=0; 22 23 FILE *fp_1,*fp_2; 24 fp_1=fopen("real_src.txt","r"); 25 if(fp_1==NULL) 26 cout<<"fopen_error"<<endl; 27 fp_2=fopen("imag_src.txt","r"); 28 if(fp_2==NULL) 29 cout<<"fopen_error"<<endl; 30 31 for(i=0;i<2048;i++) 32 { 33 fscanf(fp_1,"%lf",&u_real); 34 fscanf(fp_2,"%lf",&u_imag); 35 complex <double> u_mid(u_real,u_imag); 36 u[i]=u_mid; 37 } 38 fclose(fp_1); 39 fclose(fp_2); 40 }
从文件中读数据到变量中时，采用的函数是——fscanf()。
由于复数在定义和赋值这两个操作是同时完成的，不能分开执行，所以这里定义了一个中间变量——u_mid由它来将最终复数赋值给相应的变量——数组的元素。
1 /**************计算G=fft(g)**************/ 2 complex <double>* fft_g() 3 { 4 //得到矩阵g 5 complex <double> *g; 6 int L=czt_n+czt_m; 7 g = (complex <double> *)malloc(sizeof(complex <double>)*L); 8 for(int i=0;i<czt_n;i++) 9 g[i]=pow(czt_a,-i)*pow(czt_w,i*i/2)*u[i]; 10 for(int i=czt_n;i<L;i++) 11 g[i]=0; 12 13 //fft计算 14 fftw_complex *in_put,*out_put; 15 complex <double> *G_mid; 16 G_mid = (complex <double> *)malloc(sizeof(complex <double>)*L); 17 fftw_plan p_1; 18 int N =L; 19 in_put=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 20 out_put=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 21 for( int i=0; i < N; i++) 22 { 23 in_put[i][0] = real(g[i]); 24 in_put[i][1] = imag(g[i]); 25 //printf("%lf+%lfj ",in_put[i][0],in_put[i][1]); 26 } 27 p_1=fftw_plan_dft_1d(N,in_put,out_put, FFTW_FORWARD, FFTW_ESTIMATE); 28 fftw_execute(p_1); 29 for(int j = 0;j < N;j++) 30 { 31 complex <double> var_mid(out_put[j][0],out_put[j][1]); 32 G_mid[j]=var_mid; 33 } 34 fftw_destroy_plan(p_1); 35 fftw_free(in_put); 36 fftw_free(out_put); 37 38 return G_mid; 39 }
注意点：
1、在定义g数组时，不能使用：complex <double> g[L]这种形式，因为L是一个变量，那么数组的静态定义方式是行不通的。在动态定义中，一定要分配足够的空间，这里的g数组表示的是4096*1的矩阵空间，那么就需要malloc(sizeof(complex <double>)*L个空间。如果只是定义了一个复数的空间（complex <double> *g），那么只能存储第一个元素位置的数，其他的元素不能使用了。
2、fftw_complex类型和compelx <double>类型的数据不能强制转换，也不能隐式转换，但是可以使用实部和虚部分开来赋值：
for( int i=0; i < N; i++)
{
in_put[i][0] = real(g[i]);
in_put[i][1] = imag(g[i]);
}
for(int j = 0;j < N;j++)
{
complex <double> var_mid(out_put[j][0],out_put[j][1]);
G_mid[j]=var_mid;
}
h = zeros(L,1); %4096*1 for i=0:CZT_M-1 h(i+1,1)=CZT_W^(-i^2/2); end for i=L-CZT_N+1:L-1 h(i+1,1)=CZT_W^(-(L-i)^2/2); end H=fft(h);
注意matlab中矩阵下标是从1开始的，而在c语言中选择从0开始，要注意对应点。
从上式可以看出，h[1]——h[2048]被赋值CZT_W^(-i^2/2)，h[2050]——h[4095]被赋值为CZT_W^(-(L-i)^2/2)。h[2049]=0；
1 /**************计算H=fft(h)*********************/ 2 complex <double>* fft_h() 3 { 4 int i=0; 5 //得到矩阵h 6 complex <double> *h; 7 int L=czt_n+czt_m; 8 h = (complex <double> *)malloc(sizeof(complex <double>)*L); 9 10 //给矩阵h赋值 11 for(i=0;i<czt_m;i++) 12 h[i] = pow(czt_w,-(i*i)/2); 13 h[2048]=0; 14 for(i=czt_m+1;i<L;i++) 15 h[i]=pow(czt_w,(-(L-i)*(L-i)/2)); 16 17 //fft计算 18 fftw_complex *in_put,*out_put; 19 complex <double> *H_mid; 20 H_mid = (complex <double> *)malloc(sizeof(complex <double>)*L); 21 fftw_plan p_1; 22 int N =L; 23 in_put=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 24 out_put=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 25 for( int i=0; i < L; i++) 26 { 27 in_put[i][0] = real(h[i]); 28 in_put[i][1] = imag(h[i]); 29 //printf("%lf+%lfj ",in_put[i][0],in_put[i][1]); 30 } 31 p_1=fftw_plan_dft_1d(N,in_put,out_put, FFTW_FORWARD, FFTW_ESTIMATE); 32 fftw_execute(p_1); 33 for(int j = 0;j < N;j++) 34 { 35 complex <double> var_mid(out_put[j][0],out_put[j][1]); 36 H_mid[j]=var_mid; 37 } 38 fftw_destroy_plan(p_1); 39 fftw_free(in_put); 40 fftw_free(out_put); 41 42 return H_mid; 43 44 }
Q=H.*G;
q=ifft(Q);
1 /************计算q=ifft(Q)********************/ 2 complex <double> *ifft_Q() 3 { 4 //得到矩阵Q 5 complex <double> *H=fft_h(); //调用函数得到H矩阵 6 complex <double> *G=fft_g(); //调用函数得到G矩阵 7 int L=czt_n+czt_m; 8 complex <double> *Q; 9 Q = (complex <double> *)malloc(sizeof(complex <double>)*L); 10 for(int i=0;i<L;i++) 11 { 12 Q[i]=H[i]*G[i]; 13 } 14 15 //ifft计算 16 fftw_complex *in_put,*out_put; 17 complex <double> *q_mid; 18 q_mid = (complex <double> *)malloc(sizeof(complex <double>)*L); 19 fftw_plan p_1; 20 int N =L; 21 in_put=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 22 out_put=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 23 for( int i=0; i < L; i++) 24 { 25 in_put[i][0] = real(Q[i]); 26 in_put[i][1] = imag(Q[i]); 27 //printf("%lf+%lfj ",in_put[i][0],in_put[i][1]); 28 } 29 p_1=fftw_plan_dft_1d(N,in_put,out_put, FFTW_BACKWARD, FFTW_ESTIMATE); 30 fftw_execute(p_1); 31 for(int j = 0;j < N;j++) 32 { 33 complex <double> var_mid(out_put[j][0]/N,out_put[j][1]/N); 34 q_mid[j]=var_mid; 35 } 36 fftw_destroy_plan(p_1); 37 fftw_free(in_put); 38 fftw_free(out_put); 39 return q_mid; 40 }
值得注意的是：
p_1=fftw_plan_dft_1d(N,in_put,out_put, FFTW_BACKWARD, FFTW_ESTIMATE); FFTW_BACKWARD表示进行ifft变换，但是计算的结果没有除以点数N，所以 结果是实际ifft变换的N倍，要得到实际的ifft结果还得除以N。
所以：
for(int j = 0;j < N;j++)
{
complex <double> var_mid(out_put[j][0]/N,out_put[j][1]/N);
q_mid[j]=var_mid;
}
X=zeros(CZT_M,1); for i=0:CZT_M-1 X(i+1,1)=CZT_W^(i^2/2)*q(i+1); end y=X;
1 /***************计算y的值*******************/ 2 void y() 3 { 4 complex <double> *x,*y; 5 complex <double> *q=ifft_Q(); //调用函数得到q矩阵 6 x=(complex <double>*)malloc(sizeof(complex <double>)*czt_m); 7 y=(complex <double>*)malloc(sizeof(complex <double>)*czt_m); 8 for(int i=0;i<czt_m;i++) 9 { 10 x[i]=pow(czt_w,(i*i)/2)*q[i]; 11 } 12 y=x; 13 output_to_file(y,czt_m); 14 }
注意点：
output_to_file(y,czt_m);
函数功能是将输出矩阵——y中数据打印到一个文件中去，这样方便与matlab中结果比较，因为在终端显示器上由于点数较多，不能完全显示。这个函数需要的参数有两个：其一就是矩阵名——数组地址，其二就是矩阵长度——数组长度。完整定义如下：
1 //写一个检验函数：目的是将复数输出结果打印到文件中去，以便跟matlab结果对比 2 void output_to_file(complex <double> *name,int L) 3 { 4 complex <double> *var_name=name; 5 FILE *fp=fopen("data.txt","w+"); 6 for(int i=0;i<L;i++) 7 fprintf(fp,"%lf + %lfj ",real(var_name[i]),imag(var_name[i])); 8 fclose(fp); 9 }
所以这个函数可以检验多个中间过程的矩阵取值情况（很nice 呵呵）
最后就是主函数调用部分：
1 int main() 2 { 3 src_txt_u(); //得到原始数据点 4 complex <double> *G; 5 G=fft_g(); //得到G矩阵 6 complex <double> *H; 7 H=fft_h(); //得到H矩阵 8 y(); //H.*G,q=ifft(Q)，y=x都是在这个函数中完成的 9 return 0; 10 11 }