#include <iostream>
#include <complex>
#include <vector>
#include <cmath>
using namespace std;
// 计算傅里叶变换系数
vector<complex<double>> fft_coefficients(int N) {
vector<complex<double>> coeffs(N);
for (int k = 0; k < N; k++) {
double real_part = cos(2 * M_PI * k / N);
double imag_part = sin(-2 * M_PI * k / N);
coeffs[k] = complex<double>(real_part, imag_part);
}
return coeffs;
}
// 计算DFT
vector<complex<double>> dft(vector<double> x) {
int N = x.size();
vector<complex<double>> coeffs = fft_coefficients(N);
vector<complex<double>> X(N);
for (int k = 0; k < N; k++) {
for (int n = 0; n < N; n++) {
X[k] += x[n] * pow(coeffs[k], n);
}
}
return X;
}
// 计算希尔伯特变换
vector<double> hilbert_transform(vector<double> x) {
vector<complex<double>> X = dft(x);
int N = X.size();
vector<complex<double>> H(N);
for (int k = 0; k < N; k++) {
if (k == 0) {
H[k] = complex<double>(0, 0);
} else if (k <= N / 2) {
H[k] = complex<double>(0, -1.0 / k);
} else {
H[k] = complex<double>(0, 1.0 / (N - k));
}
}
vector<complex<double>> Y(N);
for (int k = 0; k < N; k++) {
Y[k] = X[k] * H[k];
}
vector<double> y(N);
vector<complex<double>> y_complex = dft(Y);
for (int n = 0; n < N; n++) {
y[n] = y_complex[n].imag();
}
return y;
}
int main() {
vector<double> x = {1, 2, 3, 4, 5};
vector<double> y = hilbert_transform(x);
for (int i = 0; i < y.size(); i++) {
cout << y[i] << " ";
}
cout << endl;
return 0;
}
