1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

#include "kiss_fftr.h"
#include "fftw3.h"
#include	<sndfile.h>
#include "kiss_fft.h"

void fft_d(double * in, fftw_complex* out, int nfft) {
	fftw_plan p;
	//p = fftw_plan_dft_1d(nfft, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
	// 一维实数据的DFT
	p = fftw_plan_dft_r2c_1d(nfft, in, out, FFTW_ESTIMATE);
	fftw_execute(p);
	fftw_destroy_plan(p);
}

void ifft_d(fftw_complex* in, double* out, int nfft) {
	fftw_plan p;
	//p = fftw_plan_dft_1d(nfft, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
	p = fftw_plan_dft_c2r_1d(nfft, in, out, FFTW_ESTIMATE);
	fftw_execute(p);
	fftw_destroy_plan(p);
}

void testfft() {
	int method = 1;
	//const int N = 11125;
	SF_INFO sf_info;
	SNDFILE *snd_file;
	snd_file = sf_open("p232_084_8kHz.wav", SFM_READ, &sf_info);

	float buf2[M];
	//buf2 = (double *)malloc(sf_info.frames * sizeof(double));
	//sf_readf_double(snd_file, buf2, sf_info.frames);
	printf("Sample Rate : %d\n", sf_info.samplerate);
	printf("Channels : %d\n", sf_info.channels);
	printf("Sections  : %d\n", sf_info.sections);
	printf("Frames   : %d\n", (int)sf_info.frames);
	int n_frames = (int)sf_info.frames / M;

	SNDFILE	*output_file;
	SF_INFO outsf_info;

	if (!(output_file = sf_open("sine.wav", SFM_WRITE, &sf_info)))
	{
		printf("Error : Not able to open output file.\n");
		return 1;
	}
	if (method == 0) {//方法一使用kissfft
		kiss_fft_cpx cpx_in[M];
		kiss_fft_cpx cx_out[M];
		kiss_fft_cfg cfg = kiss_fft_alloc(M, 0, NULL, NULL);


		kiss_fftr_cfg icfg = kiss_fftr_alloc(M, 1, 0, 0);
		kiss_fft_cpx freq_data[M / 2 + 1];
		kiss_fft_scalar time_data[M];
		for (int n = 0; n < n_frames; n++) {
			sf_readf_float(snd_file, buf2, M);

			for (int i = 0; i < M; i++)
			{
				cpx_in[i].r = buf2[i];
				cpx_in[i].i = 0;
			}

			kiss_fft(cfg, cpx_in, cx_out);

			/* inverse FFT */
			for (int i = 0; i < M / 2 + 1; i++)
			{
				freq_data[i].r = cx_out[i].r;
				freq_data[i].i = cx_out[i].i;
			}
			kiss_fftri(icfg, freq_data, time_data);

			float write_data[M];
			for (int i = 0; i < M; i++)
			{
				write_data[i] = time_data[i] / M;
				if (i == 0) printf("%f %f %f\n", write_data[i], buf2[i], cpx_in[i].r);
			}
			sf_write_float(output_file, write_data, M);
		}
		free(cfg);
	}
	else {//方法二：使用FFTW
		double in[M];
		fftw_complex out[M*2];

		fftw_complex freq_data[M * 2];
		double time_data[M];

		for (int n = 0; n < n_frames; n++) {
			sf_readf_float(snd_file, buf2, M);
			for (int i = 0; i < M; i++)
			{
				in[i] = buf2[i];
			}
			fft_d(in, out, M);

			/* inverse FFT */
			for (int i = 0; i < M; i++)
			{
				freq_data[i][0] = out[i][0];
				freq_data[i][1] = out[i][1];
			}
			ifft_d(freq_data, time_data, M);

			float write_data[M];
			for (int i = 0; i < M; i++)
			{
				write_data[i] = time_data[i] / M;
				if (i == 0) printf("%f %f %f\n", write_data[i], buf2[i], in[i]);
			}
			sf_write_float(output_file, write_data, M);
		}
	}
	
	sf_close(output_file);
	sf_close(snd_file);
	
}