标签：
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_w_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 11 //fft计算 12 fftw_complex *in_put,*out_put; 13 complex <double> *G_mid; 14 G_mid = (complex <double> *)malloc(sizeof(complex <double>)*2048); 15 fftw_plan p_1; 16 int N =2048; 17 in_put=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 18 out_put=(fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); 19 for( int i=0; i < N; i++) 20 { 21 in_put[i][0] = real(g[i]); 22 in_put[i][1] = imag(g[i]); 23 //printf("%lf+%lfj\n",in_put[i][0],in_put[i][1]); 24 } 25 p_1=fftw_plan_dft_1d(N,in_put,out_put, FFTW_FORWARD, FFTW_ESTIMATE); 26 fftw_execute(p_1); 27 for(int j = 0;j < N;j++) 28 { 29 complex <double> var_mid(out_put[j][0],out_put[j][1]); 30 G_mid[j]=var_mid; 31 } 32 33 for(int j = 0;j < N;j++) 34 { 35 //printf("%lf+%lfj\n",real(G_mid[j]),imag(G_mid[j])); 36 } 37 fftw_destroy_plan(p_1); 38 fftw_free(in_put); 39 fftw_free(out_put); 40 return G_mid; 41 }
注意点：
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;
}
标签：
原文地址：http://www.cnblogs.com/data1213/p/4998747.html
