DSP的一些代码实现

老师课上讲完FFT之后感叹”难道你们没有编代码的冲动吗“--遂编。

不得不说，其实还是站在前人的肩膀上，代码中最大的一关应该是复数在C++的实现，想来想去不知道该怎么写最简。后来发现有现成的complex库...感谢伟大的前人程序员。

上学期哐哐写Verilog，假期哐哐写python，重拾C++感觉熟悉又陌生。写的还是太菜了，不过记录一下吧。

#include <iostream>
#define _USE_MATH_DEFINES
#include <cmath>
#include <complex>
#include <vector>
#include <fstream>
using namespace std;
typedef complex<double> cd;
void fft(cd* arr, int N);
#define Pi    3.14159265358979323846
void DFT(cd* arr, int N);
void ifft(cd* arr, int N);
int main() {
	//DIT
	cd time[1024];
	cd freq[1024];
	int len, N;
	double re, im;



	cout << "请输入时域复序列长度" << endl;
	cin >> len;
	cout << "请输入复序列具体内容，先实部后虚部" << endl;
	for (int i = 0; i < len; i++)
	{
		cin >> re >> im;
		time[i].real(re);
		time[i].imag(im);
	}
	cout << "请输入FFT点数" << endl;
	cin >> N;

	for (int i = 0; i < len; i++)
		cout << "x[" << i << "]:" << time[i].real() << " + j" << time[i].imag() << endl;

	cout << endl;
	//不够N则填0
	for(int i =0;i<N;i++){
		if (i > len - 1)
			time[i] = 0;

		freq[i] = time[i];
	}
	//fft(freq,N);
	//DFT(freq, N);
	ifft(freq, N);
	cout << endl;

	for (int i = 0; i < N; i++){
		cout << "x["<<i<<"]:" << freq[i].real() << " + j" << freq[i].imag() << endl;

	}


	return 0;
			

};

void fft(cd * arr,int N) {


	double arg = 2 * Pi / N;
		//e的-j 2pi /N
	cd one(1, 0);
	cd temp;

	if (N != 2) {
		arg = 2 * Pi / N;

		//非2点FFT，就先奇偶分，然后N减半。
		cd* temp_a = new cd[N/2];
		cd* temp_b = new cd[N/2];

		for (int i = 0; i < N/2; i++)
		{
			temp_a[i] = arr[2 * i];
			temp_b[i] = arr[2 * i + 1];
		}

		fft(temp_a, N / 2);
		fft(temp_b, N / 2);

		for (int i = 0; i < N / 2; i++)
		{
			arg = (2 * Pi / N) * i;
			cd Wn((fabs(cos(arg)) < 1e-10)?0:cos(arg), (fabs(sin(arg)) < 1e-10) ?0: -1 * sin(arg));
			temp = temp_a[i];
			//蝶形操作
			temp_a[i] = temp_a[i] * one + Wn * temp_b[i];
			temp_b[i] = temp * one - Wn * temp_b[i];
		}
		for (int i = 0; i < N / 2; i++)
		{
			arr[i] = temp_a[i];
			arr[i + N/2] = temp_b[i];
			
		}
				
		delete[]temp_a;
		delete[]temp_b;

		return;

	}
	else {
		//2点DFT，就进行操作
		arg = (2 * Pi / N)* 0;
		cd Wn(cos(arg), -1 * sin(arg));
		temp = arr[0];
		//蝶形操作
		arr[0] = arr[0] * one + Wn * arr[1];
		arr[1] = temp * one - Wn * arr[1];

		return;
	}
		
}

void DFT(cd *arr,int N) {

	cd* tmp = new cd[N];
	double arg = 2 * Pi / N;

	for (int k = 0; k < N; k++) {
		tmp[k] = 0;
		
		for (int n = 0; n < N; n++)
		{
			cd Wn((fabs(cos(arg*k*n)) < 1e-10) ? 0 : cos(arg*k*n), (fabs(sin(arg * k * n)) < 1e-10) ? 0 : -1 * sin(arg * k*n));
			tmp[k] += Wn * arr[n];

		}

	}
	for (int i = 0; i < N; i++)
		arr[i] = tmp[i];

	delete[]tmp;
};


void ifft(cd* arr, int N) {

	//1、freq取共轭,同时/N
	for (int i = 0; i < N; i++) {
		arr[i] = conj(arr[i]);
		arr[i] /= N;
	}
	//2、fft
	fft(arr, N);

	//3、结果的时域再取共轭

	for (int i = 0; i < N; i++)
		arr[i] = conj(arr[i]);
}

结果：

IFFT/FFT互验证了一下。

（浮点那块还是有点问题，因为pi的问题，可以简单的加个if，以后再说吧）

时隔一年，自己的代码水平还是一样烂。。。菜鸡勿喷QAQ