experiment-3.cpp

计算离散信号的振幅谱和相位谱 1. 编写计算DFT的程序。 2. 编写计算振幅谱和相位谱的程序。

#include"math.h"
#include"stdio.h"
#include"stdlib.h"
void DFT(float *x, float *xr, float *xi, int N, float s);
extern void signal(float *h, float f, float dt, float m, int N);
extern void output_h(float *x, float T, int N);
extern void output_S_3(float *Xamp, float *phi, float T, int N);
extern void FFT(float *x, float *xr, float *xi, int N);
extern void F_Spectrum(float *Amp, float *phi, float *Xr, float *Xi, int N);

//-------------------------------------------------------//
//    试验-3：  计算离散信号的振幅谱和相位谱             //
//                                                       //
//  1.  编写计算DFT的程序。                              //
//                                                       //
//  2.  编写计算振幅谱和相位谱的程序。                   // 
//                                                       //
//  3.  用Excel软件绘图显示计算结果。                    //
//                                                       // 
//  4.  比较DFT和FFT的计算速度。                         //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//   编写程序的位置:                                     //
//                                                       //
//   1. 在函数void DFT()中完成DFT计算的编程              //
//                                                       //
//   2. 在主段 "计算振幅谱 Amp[i] 和 相位谱 phi[i]" 处   //
//  	完成计算振幅谱 Amp[i] 和 相位谱 phi[i]的编程     // 
//                                                       //
//-------------------------------------------------------//
void main()
{

    float x[256], h[256], Xr[256], Xi[256],
		  deltat, M, f;
	int   N;

//-------------------------------------------------------//
//             生成离散信号并计算幅频特性                //
//-------------------------------------------------------//

    deltat=float(0.002);  // 时域采样间隔(sec)
    N=256;                // 数据长度(点数)
	f=20;                 // 信号的频率
	M=float(1.1);
    signal(x, f, deltat, M, N);

//--------------------------------------------//
//       调用DFT计算离散傅立叶变换 DFT        //
//--------------------------------------------//

    DFT(x, Xr, Xi, N, 1);

//######################################################//
//------------------------------------------------------//
//  编写计算振幅谱 Amp[i] 和 相位谱 phi[i] 的程序 !!!   //
//------------------------------------------------------//
//######################################################//
    float Amp[256], phi[256];
	int i;
	for(i=0;i<N;i++)
	{
		Amp[i]=sqrt(Xr[i]*Xr[i]+Xi[i]*Xi[i]);
		phi[i]=atan(Xi[i]/Xr[i]);
	}





//------------------------------------------//
//    输出振幅谱 Amp[i] 和 相位谱 phi[i]    //
//------------------------------------------//
    output_S_3(Amp, phi, deltat, N);

//------------------------------------------//
//           检验离散傅立叶反变换           //
//------------------------------------------//
    DFT(h, Xr, Xi, N, -1);
    output_h(h, deltat, N);
    output_h(x,deltat,N);


//------------------------------------------//
//           比较DFT和FFT的计算速度         //
//------------------------------------------//

	printf("  Test for FFT&DFT !!\n\n");
	int   k;
    for(k=0; k<2000; k++)
   {    DFT(x, Xr, Xi, N, 1);  }
    {    FFT(x, Xr, Xi, N);     }


//-----------------------------------------------//
        printf("  End of caculating !!\n\n");
//-----------------------------------------------//
}



//####################################################//
//                   计算DFT的程序                    //
//----------------------------------------------------//
//         在下列函数段中编写计算DFT的程序 !!!        //
//----------------------------------------------------//
//                                                    //
//     x[]:  一维数组, 时间域的离散信号               //
//    Xr[]:  一维数组, DFT的实部                      //
//    Xi[]:  一维数组, DFT的虚部                      //
//       N:  信号长度（点数）                         //
//       s:  正反变换控制符，s=1.0 时为正变换；       //
//                           s=-1.0时为反变换         //
//                                                    //
//----------------------------------------------------//
//####################################################//
void DFT(float *x, float *Xr, float *Xi, int N, float s)
{int n,k,i,j;
 if(s=1.0)
 {for(k=0;k<N;k++)
{   Xr[k]=Xi[k]=0;
	for(n=0;n<N;n++)
	{
		Xr[k]+=x[n]*cos(2*3.14*n*k/N);
		Xi[k]+=-x[n]*sin(2*3.14*n*k/N);
	}
 }}
 if(s=-1.0)
 {for(i=0;i<N;i++)
 {x[i]=0;
    for(j=0;j<N;j++)
	{x[i]+=Xr[j]*cos(2*3.14*i*j/N)/N-Xi[j]*sin(2*3.14*i*j/N)/N;}}


 }








}