# CZT变换（chirp z-transform）

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
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 csub(comp a,comp b);
comp czero();
comp cpow(comp a,double n);
float arg(comp a);

#endif
View Code

CZT.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

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<<m)) m++;
for(i=0,j=0; i<FFTn-1; i++)
{
if(i<j)
{
t=data[j];
data[j]=data[i];
data[i]=t;
}
k=FFTn/2;
while(k<=j)
{
j-=k;
k/=2;
}
j+=k;
}
le=1;
for(l=0; l<m; l++ )
{
le*=2;
le1=le/2;
u.re=(float)1.0;
u.im=(float)0.0;
w.re=(float)cos(pi/le1);
w.im=(float)(inverse*sin(pi/le1));
for(j=0; j<le1; j++)
{
for(i=j; i<FFTn; i+=le)
{
ip=i+le1;
t.re=data[ip].re*u.re-data[ip].im*u.im;
t.im=data[ip].re*u.im+data[ip].im*u.re;
data[ip].re=data[i].re-t.re;
data[ip].im=data[i].im-t.im;
data[i].re+=t.re;
data[i].im+=t.im;
}
temp1=u.re;
temp2=u.im;
u.re=(float)(temp1*w.re-temp2*w.im);
u.im=(float)(temp1*w.im+temp2*w.re);
}
}
if(inverse>0) return;
for(i=0; i<FFTn; i++)
{
data[i].re/=FFTn;
data[i].im/=FFTn;
}
return;
}
View Code

FFT.h

#ifndef    _FFT1_H_
#define    _FFT1_H_

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

#endif
View Code

• matlab

• C

频率细化直接查找：

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 211
plot(linspace(f1,f2,length(y)),abs(y));
hold on;
scatter(fre_est,abs(y(pos)),'r*');
subplot 212
plot(t/max(t)*fs,abs(y1))

C读取txt数据：

//
#include <stdio.h>
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 <stdlib.h>
#include <stdio.h>
#include <math.h>
#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<N; i++) {
r=fscanf(f,"%f",&da);  //读取数据
x[i].re = da;
x[i].im = 0.0;
printf("n = %f\n",x[i].re );
}
fclose(f);

//

A0 = 1.0;                   //起始抽样点z0的矢量半径长度
Theta0 = 2.0*PI*fstart/fs;             //起始抽样点z0的相角
A.re = (float)(A0*cos(Theta0));
A.im = (float)(A0*sin(Theta0));

Phi0 = (fend-fstart)*2.0*PI/M/fs;      //两相邻抽样点之间的角度差
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++)
{
absx[i] = sqrt(xCZT[i].re*xCZT[i].re + xCZT[i].im*xCZT[i].im);//频谱
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);
}
printf("The Chirp-Z Transfrom square:\n");
for (i = 0 ;i<M ;i++)
{
printf("%lf\n",absx[i]);
}
//找出频谱峰值
maxcache = 0;
pos = 0;
for (i = 0 ;i<M ;i++)
{
maxdata = max_ab(maxcache,absx[i]);
if (maxdata == absx[i]) pos++;
maxcache = maxdata;
}
//频率估计
fre_est = fstart + (pos-1.0)*(fend-fstart)/M;
printf("frequency estiamton: %lf",fre_est);
//绘制频谱
//
system("pause");
}

/*
Function: max
*/
double max_ab(double x,double y)
{
return(x>y?x:y);
}

/*----------------函数说明----------------------
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]);
}

}

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));
%存入txt
fp=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];
end
end

[val,pos] = sort(dataall,'descend');
pos = posall(pos);
fre_est = (pos-1)*(f2-f1)/m+f1;
figure()

y1 = fft(sig);

subplot 211
plot(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 212
plot(t/max(t)*fs,abs(y1))
axis([95,110,0,1500]);

C：

/*=============================
Chirp-Z Transform
=============================*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#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<N; i++) {
r=fscanf(f,"%f",&da);  //读取数据
x[i].re = da;
x[i].im = 0.0;
printf("n = %f\n",x[i].re );
}
fclose(f);

//

A0 = 1.0;                   //起始抽样点z0的矢量半径长度
Theta0 = 2.0*PI*fstart/fs;             //起始抽样点z0的相角
A.re = (float)(A0*cos(Theta0));
A.im = (float)(A0*sin(Theta0));

Phi0 = (fend-fstart)*2.0*PI/M/fs;      //两相邻抽样点之间的角度差
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++)
{
absx[i] = sqrt(xCZT[i].re*xCZT[i].re + xCZT[i].im*xCZT[i].im);//频谱
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);
}
printf("The Chirp-Z Transfrom square:\n");
//for (i = 0 ;i<M ;i++)
//{
//    printf("%lf\n",absx[i]);
//}
//找出频谱峰值
//maxcache = 0;
//pos = 0;
//for (i = 0 ;i<M ;i++)
//   {
//    maxdata = max_ab(maxcache,absx[i]);
//    if (maxdata == absx[i]) pos++;
//    maxcache = maxdata;
//   }
////频率估计
//fre_est = fstart + (pos-1.0)*(fend-fstart)/M;
//printf("frequency estiamton: %lf",fre_est);
//找第二个频率分量
for (i = 0 ;i < M-1 ;i++)
{
diffdata[i] = absx[i+1]-absx[i];
}
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<flag ;i++)
{
maxdata = max_ab(maxcache,datall[i]);
if (maxdata == datall[i]) pos1 = i;
maxcache = maxdata;
}

datall[pos1] = 0;
maxcache = -1;
pos2 = 0;
for (i = 0 ;i<flag ;i++)
{
maxdata = max_ab(maxcache,datall[i]);
if (maxdata == datall[i]) pos2 = i;
maxcache = maxdata;
}
fre_est1 = fstart + (posall[pos1]-1.0)*(fend-fstart)/M;
fre_est2 = fstart + (posall[pos2]-1.0)*(fend-fstart)/M;
printf("frequency estiamton: %lf and %lf",fre_est1,fre_est2);
//绘制频谱
//
system("pause");
}

/*
Function: max
*/
double max_ab(double x,double y)
{
return(x>y?x:y);
}

/*----------------函数说明----------------------
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

posted @ 2018-05-20 13:00  桂。  阅读(15077)  评论(7编辑  收藏  举报