-
Notifications
You must be signed in to change notification settings - Fork 5
/
operator_dft_1d.cpp
80 lines (66 loc) · 2.23 KB
/
operator_dft_1d.cpp
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
// Software: Discrete Fourier Transform Opeartor 1D (use for further computation)
// Author: Hy Truong Son
// Position: PhD Student
// Institution: Department of Computer Science, The University of Chicago
// Email: [email protected], [email protected]
// Website: http://people.inf.elte.hu/hytruongson/
// Copyright 2016 (c) Hy Truong Son. All rights reserved.
#include <iostream>
#include <fstream>
#include <sstream>
#include <cstring>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <vector>
#include <set>
#include <iterator>
#include <algorithm>
#include <ctime>
#include "mex.h"
using namespace std;
void vector2matrix(double *input, int nRows, int nCols, double **output) {
for (int i = 0; i < nRows; ++i) {
for (int j = 0; j < nCols; ++j) {
output[i][j] = input[j * nRows + i];
}
}
}
void matrix2vector(double **input, int nRows, int nCols, double *output) {
for (int i = 0; i < nRows; ++i) {
for (int j = 0; j < nCols; ++j) {
output[j * nRows + i] = input[i][j];
}
}
}
void mexFunction(int nOutputs, mxArray *output_pointers[], int nInputs, const mxArray *input_pointers[]) {
if (nInputs != 1) {
std::cerr << "The number of input parameters must be exactly 1!" << std::endl;
return;
}
// Input the dimension
int N = mxGetScalar(input_pointers[0]);
// Memory allocation
double **Re_w = new double* [N];
double **Im_w = new double* [N];
for (int u = 0; u < N; ++u) {
Re_w[u] = new double [N];
Im_w[u] = new double [N];
}
// Computation
for (int u = 0; u < N; ++u) {
double t = 2.0 * M_PI * (double)(u) / (double)(N);
for (int v = u; v < N; ++v) {
double alpha = t * (double)(v);
Re_w[u][v] = cos(alpha);
Im_w[u][v] = - sin(alpha);
Re_w[v][u] = Re_w[u][v];
Im_w[v][u] = Im_w[u][v];
}
}
// Return results
output_pointers[0] = mxCreateDoubleMatrix(N, N, mxREAL);
matrix2vector(Re_w, N, N, mxGetPr(output_pointers[0]));
output_pointers[1] = mxCreateDoubleMatrix(N, N, mxREAL);
matrix2vector(Im_w, N, N, mxGetPr(output_pointers[1]));
}