CZT变换（chirp z-transform）

作者：桂。

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

链接：http://www.cnblogs.com/xingshansi/p/9063131.html 

---

前言

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

一、代码实现

原理主要参考MATLAB接口：

对应C代码实现：

Complex.c

/*=============================
                          Chirp-Z Transform
=============================*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#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<N; i++)
    {
        printf("%10.4f",x[i].re);
        printf("%10.4f\n",x[i].im);
    }

    printf("The Chirp-Z Transfrom:\n");
    for (i = 0 ;i<M ;i++)
    {
        printf("%10.4f",xCZT[i].re);
        printf("%10.4f\n",xCZT[i].im);
    }

}


/*----------------函数说明----------------------
Name: CZT
Function: Chirp-Z Transform
Para:  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<N+M-1);
    L = i;

    h = (comp*)calloc(L,sizeof(comp));
    g = (comp*)calloc(L,sizeof(comp));
    pComp = (comp*)calloc(L,sizeof(comp));

    for (i = 0; i<N; i++)
    {
        tmp1 = cpow(A,-i);
        tmp2 = cpow(W, i*i/2.0);
        tmp = cmul(tmp1,tmp2);
        g[i] = cmul(tmp,x[i]);
    }
    for (i = N;i<L; i++)
    {
        g[i] =czero();
    }

    FFT(g,L,1);

    for (i = 0;i<=M-1;i++)
    {
        h[i] = cpow(W, -i*i/2.0);
    }
    for (i=M; i<=L-N;i++)
    {
        h[i] =czero();
    }
    for (i = L-N+1; i<=L;i++)
    {
        h[i] = cpow(W,-(L-i)*(L-i)/2.0);
    }

    FFT(h,L,1);

    for (i = 0; i<L; i++)
    {
        pComp[i] = cmul(h[i],g[i]);
    }

    FFT(pComp,L,-1);         //IDFT

    for (i = 0; i<M;i++)
    {
        tmp = cpow(W,i*i/2.0);
        xCZT[i] = cmul(tmp,pComp[i]);
    }


}

View Code

Complex.h

/*===========================
Define comp as complex type
cmplx     c = (a,b)
cmul     c=a*b
conjg    c=a'
cabs1    f=|a|
cabs2    f=|a|**2
cadd     c=a+b
csub     c=a-b
czero    c=(0.0,0.0)
===========================*/
#ifndef COMPLEX_H
#define COMPLEX_H
#include <math.h>
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

View Code

CZT.c

/*=============================