作者：桂。

时间：2018-05-20  12:04:24

链接： 

---

前言

相比DFT，CZT是完成频谱细化的一种思路，本文主要记录CZT的C代码实现。

一、代码实现

原理主要参考MATLAB接口：

对应C代码实现：

Complex.c

/*=============================                          Chirp-Z Transform=============================*/#include 
      
       #include 
       
        #include 
        
         #include "Complex.h"#include "FFT.h"void CZT(comp* x, int  N, comp A, comp W, comp *xCZT, int  M);void main(){    int i;    int N,M;    double PI;    double A0,Theta0;    double W0,Phi0;    comp* x;    comp* xCZT;    comp A,W;    PI = 3.1415926;    N = 5;              //信号长度    M = 10;            //chirp-z 变换输出长度    x = (comp *)calloc(N,sizeof(comp));    xCZT  = (comp *)calloc(M,sizeof(comp));    for (i = 0;i < N; i++)    {        x[i].re=(float)(i-3);        x[i].im = 0.0;    }    A0 = 1.0;                   //起始抽样点z0的矢量半径长度    Theta0 = 0.0;             //起始抽样点z0的相角    A.re = (float)(A0*cos(Theta0));    A.im = (float)(A0*sin(Theta0));    Phi0 = 2.0*PI/M;      //两相邻抽样点之间的角度差    W0 = 1.0;                  //螺线的伸展率    W.re = (float)(W0*cos(-Phi0));    W.im = (float)(W0*sin(-Phi0));    CZT(x,N,A,W,xCZT,M);    printf("The Original Signal:\n");    for (i = 0; i
        
       
      

Complex.h

/*===========================Define comp as complex typecmplx     c = (a,b)cmul     c=a*bconjg    c=a'cabs1    f=|a|cabs2    f=|a|**2cadd     c=a+bcsub     c=a-bczero    c=(0.0,0.0)===========================*/#ifndef COMPLEX_H#define COMPLEX_H#include 
      
       typedef  struct xy{    float re;    float im;}comp;comp cmplx(float a,float b);comp cmul(comp a,comp b);comp conjg(comp a);float cabs1(comp a);float cabs2(comp a);comp cadd(comp a,comp b);comp csub(comp a,comp b);comp czero();comp cpow(comp a,double n);float arg(comp a);#endif
      

CZT.c

/*=============================                          Chirp-Z Transform=============================*/#include 
      
       #include 
       
        #include 
        
         #include "Complex.h"#include "FFT.h"void CZT(comp* x, int  N, comp A, comp W, comp *xCZT, int  M);void main(){    int i;    int N,M;    double PI;    double A0,Theta0;    double W0,Phi0;    comp* x;    comp* xCZT;    comp A,W;    PI = 3.1415926;    N = 5;              //信号长度    M = 10;            //chirp-z 变换输出长度    x = (comp *)calloc(N,sizeof(comp));    xCZT  = (comp *)calloc(M,sizeof(comp));    for (i = 0;i < N; i++)    {        x[i].re=(float)(i-3);        x[i].im = 0.0;    }    A0 = 1.0;                   //起始抽样点z0的矢量半径长度    Theta0 = 0.0;             //起始抽样点z0的相角    A.re = (float)(A0*cos(Theta0));    A.im = (float)(A0*sin(Theta0));    Phi0 = 2.0*PI/M;      //两相邻抽样点之间的角度差    W0 = 1.0;                  //螺线的伸展率    W.re = (float)(W0*cos(-Phi0));    W.im = (float)(W0*sin(-Phi0));    CZT(x,N,A,W,xCZT,M);    printf("The Original Signal:\n");    for (i = 0; i
        
       
      

FFT.c

#include "FFT.h"void FFT(comp *data,int FFTn,int inverse){    comp      u,w,t;    double     temp1,temp2;    double     pi;    int       i,j,k,l,ip;    int       le,le1,m;    pi=3.1415926f;    m=0;    while(FFTn != (0x0001<
      
       0) return;    for(i=0; i
      

FFT.h

#ifndef    _FFT1_H_#define    _FFT1_H_#include "Complex.h"//========================================//功能:    实现FFT//输入:    data[in][out]: 数据指针;  FFTn[in]:FFT点数; //            inverse[in]:正反FFT标志位:1,正FFT;-1,逆FFTvoid FFT(comp *data,int FFTn,int inverse);#endif

二、仿真测试

参数设置：

结果对比：

• matlab

• C

 

利用CZT对信号进行频率估计：

 频率细化直接查找：

仿真code:

clc;clear all;close all;fs = 1000;f0 = 201.3;t = [0:199]/fs;sig = (sin(2*pi*t*f0));% %存入txt% fp=fopen('data.txt','a');%'A.txt'为文件名；'a'为打开方式：在打开的文件末端添加数据，若文件不存在则创建。% fprintf(fp,'%f ',sig);%fp为文件句柄，指定要写入数据的文件。注意：%d后有空格。% fclose(fp);%关闭文件。f1 = 190;f2 = 210;m = 100;w = exp(-1j*2*pi*(f2-f1)/(m*fs));a = exp(1j*2*pi*f1/fs);y = czt(sig,m,w,a);[val,pos] = max(abs(y));fre_est = (pos-1)*(f2-f1)/m+f1;y1 = fft(sig);figure()subplot 211plot(linspace(f1,f2,length(y)),abs(y));hold on;scatter(fre_est,abs(y(pos)),'r*');subplot 212plot(t/max(t)*fs,abs(y1))

　　C读取txt数据：

//#include 
     
        int n,r;  double d;  FILE *f;  void main() {      f=fopen("data.txt","r");      n=0;      while (1) {          r=fscanf(f,"%lf",&d);          if (1==r) {              n++;              printf("[%d]==%lg\n",n,d);          } else if (0==r) {              fscanf(f,"%*c");          } else break;      }      fclose(f);      system("pause");}
     

C语言仿真：

/*=============================                          Chirp-Z Transform=============================*/#include 
     
      #include 
      
       #include 
       
        #include "Complex.h"#include "FFT.h"#define Size 100 //信号长度void CZT(comp* x, int  N, comp A, comp W, comp *xCZT, int  M);double max_ab(double x,double y);//求最大值void main(){    int i;    int r;    int N,M;    float da;    double fre_est;    FILE *f;     double PI;    double A0,Theta0;    double W0,Phi0;    int pos;    double maxdata;    double maxcache;    float fs;    float fstart;    float fend;    comp* x;    double absx[Size];//存储绝对值    comp* xCZT;    comp A,W;    fs = 1000.0;//frequency sample    fstart = 200.0;    fend = 215.0;    PI = 3.1415926;    f=fopen("data.txt","r");      N = Size;              //信号长度    M = Size;            //chirp-z 变换输出长度    x = (comp *)calloc(N,sizeof(comp));    xCZT  = (comp *)calloc(M,sizeof(comp));     for (i = 0; i
        
         y?x:y); } /*----------------函数说明----------------------Name: CZTFunction: Chirp-Z TransformPara:  x[in][out]:待变换信号                            N[in]:信号长度         A[in]:                                                       W[in]:         M[in]:Chirp-Z变换输出长度--------------------------------------------*/void CZT(comp* x, int  N, comp A, comp W, comp* xCZT, int  M){    int i;    int L;    comp* h;    comp* g;    comp* pComp;    comp tmp,tmp1,tmp2;    i=1;    do     {        i*=2;    } while (i
        
       
      
     

原始数据：

估计结果：

与MATLAB结果一致：

如果存在两个信号呢？

MATLAB

clc;clear all;close all;fs = 256;% f0 = 202;% f1 = 208;t = [0:625]/fs;sig = single(2*cos(2*pi*102*t)+5*cos(2*pi*105*t));%存入txtfp=fopen('data.txt','a');%'A.txt'为文件名；'a'为打开方式：在打开的文件末端添加数据，若文件不存在则创建。fprintf(fp,'%f\n',sig);%fp为文件句柄，指定要写入数据的文件。注意：%d后有空格。fclose(fp);%关闭文件。f1 = 100;f2 = 108;m = 625;w = exp(-1j*2*pi*(f2-f1)/(m*fs));a = exp(1j*2*pi*f1/fs);y = czt(sig,m,w,a);data_cache = diff(abs(y));vals = abs(y);dataall = [];posall = [];flag = 0;for i = 1:length(y)-2    if ((data_cache(1,i)*data_cache(1,i+1))<0)        flag = flag + 1;        dataall = [ dataall vals(i)];        posall = [posall,i];    endend[val,pos] = sort(dataall,'descend');pos = posall(pos);fre_est = (pos-1)*(f2-f1)/m+f1;figure()y1 = fft(sig);subplot 211plot(linspace(f1,f2,length(y)),abs(y));axis([95,110,0,1500]);hold on;scatter(fre_est(1:2),ones(1,2),'r*');hold on;% scatter(fre_est,abs(y(pos)),'r*');subplot 212plot(t/max(t)*fs,abs(y1))axis([95,110,0,1500]);

　　C：

运行中如果报错：

 

需要将fopen的data.txt添加绝对路径，注意：\\而不是\。

/*=============================                          Chirp-Z Transform=============================*/#include 
      
       #include 
       
        #include 
        
         #include "Complex.h"#include "FFT.h"#define Size 626 //信号长度void CZT(comp* x, int  N, comp A, comp W, comp *xCZT, int  M);double max_ab(double x,double y);//求最大值void main(){    int i;    int r;    int flag = 0;    int N,M;    float da;    double fre_est1, fre_est2;    FILE *f;     double datall[Size] = {
    0};//新建足够大的数据    int posall[Size] = {
    0};    double diffdata [Size-1];    double PI;    double A0,Theta0;    double W0,Phi0;    int pos1,pos2;    double maxdata;    double maxcache;    float fs;    float fstart;    float fend;    comp* x;    double absx[Size];//存储绝对值    comp* xCZT;    comp A,W;    fs = 256.0;//frequency sample    fstart = 100.0;    fend = 108.0;    PI = 3.1415926;    f=fopen("data.txt","r");      N = Size;              //信号长度    M = Size;            //chirp-z 变换输出长度    x = (comp *)calloc(N,sizeof(comp));    xCZT  = (comp *)calloc(M,sizeof(comp));     for (i = 0; i
         
          < M-2 ;i++)    {        if((diffdata[i]*diffdata[i+1])<0.0)        {            datall[flag] = absx[i];            posall[flag] = i;            flag ++;        }    }    maxcache = -1;    pos1 = 0;    for (i = 0 ;i
          
           y?x:y); } /*----------------函数说明----------------------Name: CZTFunction: Chirp-Z TransformPara:  x[in][out]:待变换信号                            N[in]:信号长度         A[in]:                                                       W[in]:         M[in]:Chirp-Z变换输出长度--------------------------------------------*/void CZT(comp* x, int  N, comp A, comp W, comp* xCZT, int  M){    int i;    int L;    comp* h;    comp* g;    comp* pComp;    comp tmp,tmp1,tmp2;    i=1;    do     {        i*=2;    } while (i
          
         
        
       
      

另外：quantusII有对应的DSP builder，Vivado/ISE是否有类似工具，直接代码转HDL？ HLS算是一个，其他呢？有时间可以将*.c转成*.v试验。

一个基本应用，频率精确测量：