FFTW3Backend.cpp 4.29 KB
Newer Older
Jonas Stienen's avatar
Jonas Stienen committed
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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
/*
 *  ITAFFT, eine Wrapper-Bibliothek für schnelle Fouriertransformationen
 *
 *  Autor: Frank Wefers (Frank.Wefers@akustik.rwth-aachen.de)
 *
 *  (c) Copyright Institut für Technische Akustik (ITA), RWTH Aachen
 *
 */

// $Id: FFTW3Backend.cpp,v 1.3 2010-01-18 17:19:01 fwefers Exp $

#include "FFTW3Backend.h"

#include <cassert>
#include <ITACriticalSection.h>
#include <ITAFFT.h>
#include <ITAException.h>

ITACriticalSection csFFTW3BackendPlannerLock;

FFTW3Realization::FFTW3Realization(int type, int size, float* in, float* out, unsigned int uiFlags)
: m_plan(NULL), m_type(type), m_size(size), m_out(out)
{
	if ((size <= 0) || (in == NULL) || (out == NULL)) ITA_EXCEPT0(INVALID_PARAMETER);

	fftw_iodim dim, howmany_dim;

	csFFTW3BackendPlannerLock.enter();
	
	// [fwe 2010-07-20] Wichtiges Bugfix. Für out-of-place Transformationen die Eingabe unberührt lassen.
	m_inplace = (in == out);
	if (m_inplace) uiFlags |= FFTW_PRESERVE_INPUT;

	switch (type) {
	case ITAFFT::FFT_R2C:
		m_plan = fftwf_plan_dft_r2c_1d(size, in, (fftwf_complex*) out, uiFlags);
		m_sInfo = "FFT_R2C [fftw3]";
		break;

	case ITAFFT::FFT_C2C:
		m_plan = fftwf_plan_dft_1d(size, (fftwf_complex*) in, (fftwf_complex*) out, FFTW_FORWARD, uiFlags);
		m_sInfo = "FFT_C2C [fftw3]";
		break;

	case ITAFFT::IFFT_C2R:
		m_plan = fftwf_plan_dft_c2r_1d(size, (fftwf_complex*) in, out, uiFlags);
		m_sInfo = "IFFT_C2R [fftw3]";
		break;

	case ITAFFT::IFFT_C2C:
		m_plan = fftwf_plan_dft_1d(size, (fftwf_complex*) in, (fftwf_complex*) out, FFTW_BACKWARD, uiFlags);
		m_sInfo = "IFFT_C2C [fftw3]";
		break;

    case ITAFFT::SPLIT_FFT_R2C:
		dim.n = size;
		dim.is = dim.os = 1;

		howmany_dim.n = 1;
		howmany_dim.is = howmany_dim.os = 1;

		m_plan = fftwf_plan_guru_split_dft_r2c(1, &dim, 1, &howmany_dim, in, out, out+(size/2)+1, uiFlags);
		m_sInfo = "SPLIT_FFT_R2C [fftw3]";
		break;

	case ITAFFT::SPLIT_IFFT_C2R:
		dim.n = size;
		dim.is = dim.os = 1;

		howmany_dim.n = 1;
		howmany_dim.is = howmany_dim.os = 1;

		m_plan = fftwf_plan_guru_split_dft_c2r(1, &dim, 1, &howmany_dim, in, in+(size/2)+1, out, uiFlags);
		m_sInfo = "SPLIT_IFFT_C2R [fftw3]";
		break;

	case ITAFFT::NORMALIZED_SPLIT_IFFT_C2R:
		dim.n = size;
		dim.is = dim.os = 1;

		howmany_dim.n = 1;
		howmany_dim.is = howmany_dim.os = 1;

		m_plan = fftwf_plan_guru_split_dft_c2r(1, &dim, 1, &howmany_dim, in, in+(size/2)+1, out, uiFlags);
		m_sInfo = "NORMALIZED_SPLIT_IFFT_C2R [fftw3]";
		break;
	}

	csFFTW3BackendPlannerLock.leave();

	if (m_plan == NULL) ITA_EXCEPT0(UNKNOWN);
}

FFTW3Realization::~FFTW3Realization() {
	if (m_plan != NULL) fftwf_destroy_plan(m_plan);
}

void FFTW3Realization::execute() {
	fftwf_execute(m_plan);

	// Sonderfall NORMALIZED_SPLIT_IFFT_C2R: Hier noch die Normalisierung durchführen
	if (m_type == ITAFFT::NORMALIZED_SPLIT_IFFT_C2R)
		for (int i=0; i<m_size; i++) m_out[i] /= (float) m_size;
}

void FFTW3Realization::execute(float* in, float* out) {
	if (m_inplace) {
		if (in != out) ITA_EXCEPT1(INVALID_PARAMETER, "Attempt to execute in-place transform out-of-place");
	} else {
		if (in == out) ITA_EXCEPT1(INVALID_PARAMETER, "Attempt to execute out-of-place transform in-place");
	}

	switch (m_type) {
	case ITAFFT::FFT_R2C:
		fftwf_execute_dft_r2c(m_plan, in, (fftwf_complex*) out);
		break;

	case ITAFFT::IFFT_C2R:
		fftwf_execute_dft_c2r(m_plan, (fftwf_complex*) in, out);
		break;

	case ITAFFT::FFT_C2C:
	case ITAFFT::IFFT_C2C:
		fftwf_execute_dft(m_plan, (fftwf_complex*) in, (fftwf_complex*) out);
		break;

	case ITAFFT::SPLIT_FFT_R2C:
		fftwf_execute_split_dft_r2c(m_plan, in, out, (out+(m_size/2)+1));
		break;

	case ITAFFT::SPLIT_IFFT_C2R:
		fftwf_execute_split_dft_c2r(m_plan, in, (in+(m_size/2)+1), out);
		break;

	case ITAFFT::NORMALIZED_SPLIT_IFFT_C2R:
		fftwf_execute_split_dft_c2r(m_plan, in, (in+(m_size/2)+1), out);
		
		// Normalisierung durchführen
		for (int i=0; i<m_size; i++) out[i] /= (float) m_size;
		
		break;
	}
}

std::string FFTW3Realization::toString() {
	return m_sInfo;
}

// -----------------------------------

FFTW3Backend* FFTW3Backend::m_pInstance = NULL;

FFTW3Backend* FFTW3Backend::getInstance() {
	if (m_pInstance == NULL) m_pInstance = new FFTW3Backend();
	return m_pInstance;
}

ITAFFTRealization* FFTW3Backend::plan(int type, int size, float* in, float* out, unsigned int uiFlags) {
	return new FFTW3Realization(type, size, in, out, uiFlags);
}