标签：https 技术分享 说明 关闭 call opened size abs text
作者：桂。
时间：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]); } }
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
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]); } }
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; }
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
二、仿真测试
参数设置：
结果对比：
利用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 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))
一个基本应用，频率精确测量：一种卫星信号载波频率精确估计算法.pdf
标签：https 技术分享 说明 关闭 call opened size abs text
原文地址：https://www.cnblogs.com/xingshansi/p/9063131.html
