滤波.cpp

针对基于DSP的使用C++编译的数字信号的高低频滤波的编程源代码

#include <iostream>
#include <fstream>
#include <math.h>
#include <iomanip>
#include "nrutil.h"

#define   N   256

using namespace std;

//傅立叶变换函数
void Gfft(double *ax, double *ay, int count,  int fftType)
{
	
    int mm=int(log(count)/log(2)+0.9999); 
    count=(int)pow(2,mm);
    int n = count;
    int n2 = n;
    for (int k = 1;k<mm+1;k++)
	{
		int n1 = n2;
		n2 = n2 / 2;
		double  e = 6.28318530717959f*fftType / n1;
		double  a1 = 0.0;
		for (int j = 1;j<n2+1;j++)
		{
			double cc = (double)cos(a1);
			double ss = (double)sin(a1);
			a1 = j * e;
			for (int i = j;i<n+1;i+=n1)
			{
				int q = i + n2;
				double xt = ax[i] - ax[q];
				ax[i] = ax[i] + ax[q];
				double  yt = ay[i] - ay[q];
				ay[i] = ay[i] + ay[q];
				ax[q] = cc * xt + ss * yt;
				ay[q] = cc * yt - ss * xt;
			}
		}
    }
    int j = 1;
    int n1 = n - 1;
    for(int i=1;i<n1+1;i++)
	{
		if (i < j )
		{
			double xt = ax[j];
			ax[j] = ax[i];
			ax[i] = xt;
			xt = ay[j];
			ay[j] = ay[i]; 
			ay[i] = xt;
		}
		k = n / 2;
		while (k < j)
		{
			j = j - k;
			k = k / 2;
		}
		j = j + k;
		if(fftType==-1)
		{
			ax[i]=ax[i]/count;
			ay[i]=ay[i]/count;
		}
    }
}

void main()
{
	int           i;
	double        tmp=0.0;
	double        *A,*ax,*ay,*az,*aa;

	//读取原数据
	A=dvector(1,N);
	aa=dvector(1,N);
	ax=dvector(1,N);
	ay=dvector(1,N);
	az=dvector(1,N);
	fstream       infile("Data.dat",ios_base::in);
	if (!infile)
	{
		cout<<"不能打开 Data.dat ！"<<endl;
		abort();
	}
	for (i=1;i<=N;i++)
	{
		infile>>tmp;
		A[i]=tmp;
	}
	infile.close();

	for (i=1;i<=N;i++)
	{
		ax[i]=A[i];
		ay[i]=0.0;
	}
	Gfft(ax,ay,N,1);
	for (i=1;i<=N;i++)
	{
//		az[i]=sqrt(ay[i]*ay[i]+ax[i]*ax[i]);   //输出的原信号的频谱
	}

	/*低通滤波
	for (i=19;i<=N;i++)
	{
		ax[i]=0.0;
		ay[i]=0.0;
	}
*/	Gfft(ax,ay,N,-1);

	Gfft(ax,ay,N,1);
	for (i=1;i<=N;i++)
	{
		az[i]=sqrt(ay[i]*ay[i]+ax[i]*ax[i]);
	}

/*	
	//时域滤波
	for (i=1;i<=20;i++)
	{
		aa[i]=1;
		az[i]=0.0;
	}
	for (;i<=N;i++)
	{
		aa[i]=0.0;
		az[i]=0.0;
	}
	Gfft(aa,az,N,-1);
    aa[N]=0.0;
	double tp[N+1],tpp[N+1];

	for (i=1;i<=N;i++)
	{
		tp[i]=0.0;
		for (int j=i;j<=N;j++)
		{
			tp[i]+=aa[j]*A[N-j+i];
			//tpp[i]+=az[j]*ay[N-j+i];
		}
	}
	
*/	

	/*高通滤波
	for (i=180;i<=N;i++)
	{
		ax[i]=0.0;
		ay[i]=0.0;
	}
	
	Gfft(ax,ay,N,-1);
*/	Gfft(ax,ay,N,1);
	for (i=1;i<=N;i++)
	{
		az[i]=sqrt(ay[i]*ay[i]+ax[i]*ax[i]);
	}

	
	//带通滤波
	for (i=1;i<=100;i++)
	{
		ax[i]=0.0;
		ay[i]=0.0;
	}
	for (i=160;i<=N;i++)
	{
		ax[i]=0.0;
		ay[i]=0.0;
	}
	
	Gfft(ax,ay,N,-1);
	Gfft(ax,ay,N,1);
	for (i=1;i<=N;i++)
	{
		az[i]=sqrt(ay[i]*ay[i]+ax[i]*ax[i]);
	}


	//输出数据
	fstream       outfile("Out.dat",ios_base::out);
	if (!outfile)
	{
		cout<<"不能打开 Out.dat ！"<<endl;
		abort();
	}
	for (i=1;i<N;i++)
	{
		outfile<<setiosflags(ios_base::fixed)
		       <<setprecision(9)<<ax[i]<<endl;
	}
	outfile.close();

	//数组释放内存空间
	free_dvector(A,1,N);
	free_dvector(aa,1,N);
	free_dvector(ax,1,N);
	free_dvector(ay,1,N);
	free_dvector(az,1,N);

	cout<<endl<<"\t程序运行成功！"<<endl<<endl;
}