代码拉取完成，页面将自动刷新
#pragma once
#include <iostream>
#include <vector>
#include <complex>
#include <functional>
#include <omp.h>
namespace FFT {
using namespace std;
typedef complex<double> cpx;
typedef vector<double> vector1d;
typedef vector<vector<double>> vector2d;
typedef vector<vector<vector<double>>> vector3d;
typedef vector<cpx> cpx1d;
typedef vector<vector<cpx>> cpx2d;
typedef vector<vector<vector<cpx>>> cpx3d;
const double PI = acos(-1);
inline int __expand_exp2(int n) { return (int)pow(2, ceil(log2((double)n))); } //扩充到2的指数
inline bool __is_exp2(int n) { return pow(2, static_cast<int>(log2((double)n))) == n; }
inline int __get_reverse_bin(int num, int n) { //获得反比特序的数字
int ret = 0;
for (int i = 0; i < n; i++)
{
if (num & (1 << i)) ret |= (1 << (n - 1 - i));
}
return ret;
}
/*初始化数组*/
template<typename T> void init(vector<T>& vec1d, int dim1, T val);
template<typename T> void init(vector<vector<T>>& vec2d, int dim1, int dim2, T val); // 行x列
template<typename T> void init(vector<vector<vector<T>>>& vec3d, int dim1, int dim2, int dim3, T val); //行x列x深度
/*把大小扩充到2的指数大小，如：63->64 = 2 ^ 6； 11->16 = 2 ^ 4。多余的部分用val填充*/
template<typename T> void expand(vector<T>& vec1d, T val);
template<typename T> void expand(vector<vector<T>>& vec2d, T val);
template<typename T> void expand(vector<vector<vector<T>>>& vec3d, T val);
/*打印数组*/
void print(vector1d& vec1d);
void print(vector2d& vec2d);
void print(vector3d& vec3d);
void print(cpx1d& vec1d);
void print(cpx2d& vec2d);
void print(cpx3d& vec3d);
/*实数和复数转换*/
void real2cpx(vector1d& src, cpx1d& dst);
void real2cpx(vector2d& src, cpx2d& dst);
void real2cpx(vector3d& src, cpx3d& dst);
void cpx2real(cpx1d& src, vector1d& dst);
void cpx2real(cpx2d& src, vector2d& dst);
void cpx2real(cpx3d& src, vector3d& dst);
/*对每一个元素的操作*/
template<typename T> void dosome(vector<T>& vec1d, function<void(T&)> func);
template<typename T> void dosome(vector<vector<T>>& vec2d, function<void(T&)> func);
template<typename T> void dosome(vector<vector<vector<T>>>& vec3d, function<void(T&)> func);
/*傅里叶变换*/
void fft(cpx1d& src, cpx1d& dst, bool inverse = false, double freq = 2 * FFT::PI); //inverse: false:逆傅里叶变换 true:正傅里叶变换
void fft2d(cpx2d& src, cpx2d& dst, bool inverse = false, double freq = 2 * FFT::PI); //inverse: false:逆傅里叶变换 true:正傅里叶变换
void fft3d(cpx3d& src, cpx3d& dst, bool inverse = false, double freq = 2 * FFT::PI); //inverse: false:逆傅里叶变换 true:正傅里叶变换
/*转换成2维频谱图*/
void spectrum(cpx2d& src, vector2d& dist);
}
template<typename T>
void FFT::init(vector<T>& vec1d, int dim1, T val)
{
vec1d.resize(dim1);
fill(vec1d.begin(), vec1d.end(), val);
}
template<typename T>
void FFT::init(vector<vector<T>>& vec2d, int dim1, int dim2, T val)
{
vec2d.resize(dim1);
#pragma omp parallel for
for (int i = 0; i < dim1; ++i) {
vec2d[i].resize(dim2);
fill(vec2d[i].begin(), vec2d[i].end(), val);
}
}
template<typename T>
void FFT::init(vector<vector<vector<T>>>& vec3d, int dim1, int dim2, int dim3, T val)
{
vec3d.resize(dim1);
#pragma omp parallel for
for (int i = 0; i < dim1; ++i) {
vec3d[i].resize(dim2);
for (int j = 0; j < dim2; ++j) {
vec3d[i][j].resize(dim3);
fill(vec3d[i][j].begin(), vec3d[i][j].end(), val);
}
}
}
void FFT::print(vector1d& vec1d)
{
size_t dim1 = vec1d.size();
cout << "[";
for (int i = 0; i < dim1; ++i) {
if (i == dim1 - 1) cout << vec1d[i];
else cout << vec1d[i] << ", ";
}
cout << "]" << endl;
}
void FFT::print(vector2d& vec2d)
{
size_t dim1 = vec2d.size();
cout << "[";
for (int i = 0; i < dim1; ++i) {
print(vec2d[i]);
}
cout << "]" << endl;
}
void FFT::print(vector3d& vec3d)
{
size_t dim1 = vec3d.size();
cout << "[";
for (int i = 0; i < dim1; ++i) {
print(vec3d[i]);
}
cout << "]" << endl;
}
void FFT::print(cpx1d& vec1d)
{
size_t dim1 = vec1d.size();
cout << "[";
for (int i = 0; i < dim1; ++i) {
if (i == dim1 - 1) cout << "(" << vec1d[i].real() << "+" << vec1d[i].imag() << "j)";
else cout << "(" << vec1d[i].real() << "+" << vec1d[i].imag() << "j), ";
}
cout << "]" << endl;
}
void FFT::print(cpx2d& vec2d)
{
size_t dim1 = vec2d.size();
cout << "[";
for (int i = 0; i < dim1; ++i) {
print(vec2d[i]);
}
cout << "]" << endl;
}
void FFT::print(cpx3d& vec3d)
{
size_t dim1 = vec3d.size();
cout << "[";
for (int i = 0; i < dim1; ++i) {
print(vec3d[i]);
}
cout << "]" << endl;
}
void FFT::real2cpx(vector1d& src, cpx1d& dst)
{
size_t dim1 = src.size();
dst.resize(dim1);
#pragma omp parallel for
for (int i = 0; i < dim1; ++i) {
dst[i] = cpx(src[i], 0.0);
}
}
void FFT::real2cpx(vector2d& src, cpx2d& dst)
{
size_t dim1 = src.size(), dim2 = src[0].size();
dst.resize(dim1);
#pragma omp parallel for
for (int i = 0; i < dim1; ++i) {
dst[i].resize(dim2);
for (int j = 0; j < dim2; ++j) {
dst[i][j] = cpx(src[i][j], 0.0);
}
}
}
void FFT::real2cpx(vector3d& src, cpx3d& dst)
{
size_t dim1 = src.size(), dim2 = src[0].size(), dim3 = src[0][0].size();
dst.resize(dim1);
#pragma omp parallel for
for (int i = 0; i < dim1; ++i) {
dst[i].resize(dim2);
for (int j = 0; j < dim2; ++j) {
dst[i][j].resize(dim3);
for (int k = 0; k < dim3; ++k) {
dst[i][j][k] = cpx(src[i][j][k], 0.0);
}
}
}
}
void FFT::cpx2real(cpx1d& src, vector1d& dst)
{
size_t dim1 = src.size();
dst.resize(dim1);
#pragma omp parallel for
for (int i = 0; i < dim1; ++i) {
dst[i] = sqrt(src[i].real() * src[i].real() + src[i].imag() * src[i].imag());
}
}
void FFT::cpx2real(cpx2d& src, vector2d& dst)
{
size_t dim1 = src.size(), dim2 = src[0].size();
dst.resize(dim1);
#pragma omp parallel for
for (int i = 0; i < dim1; ++i) {
dst[i].resize(dim2);
for (int j = 0; j < dim2; ++j) {
dst[i][j] = sqrt(src[i][j].real() * src[i][j].real() + src[i][j].imag() * src[i][j].imag());
}
}
}
void FFT::cpx2real(cpx3d& src, vector3d& dst)
{
size_t dim1 = src.size(), dim2 = src[0].size(), dim3 = src[0][0].size();
dst.resize(dim1);
#pragma omp parallel for
for (int i = 0; i < dim1; ++i) {
dst[i].resize(dim2);
for (int j = 0; j < dim2; ++j) {
dst[i][j].resize(dim3);
for (int k = 0; k < dim3; ++k) {
dst[i][j][k] = sqrt(src[i][j][k].real() * src[i][j][k].real() + src[i][j][k].imag() * src[i][j][k].imag());
}
}
}
}
void FFT::fft(cpx1d& src, cpx1d& dst, bool inverse, double freq)
{
if (!__is_exp2((int)src.size())) {
FFT::expand(src, FFT::cpx(0, 0));
}
int lim = (int)src.size();
int opt = inverse ? -1 : 1;
// https://blog.csdn.net/Isaac320/article/details/122255992
FFT::cpx1d tempA(lim);
for (int i = 0; i < lim; i++)
{
tempA[i] = src[FFT::__get_reverse_bin(i, (int)log2(lim))];
}
FFT::cpx1d WN(lim / 2);
//生成WN表,避免重复计算
for (int i = 0; i < lim / 2; i++)
{
WN[i] = FFT::cpx(cos(freq * i / lim), opt * -sin(freq * i / lim));
}
//蝶形运算
for (int steplenght = 2; steplenght <= lim; steplenght *= 2)
{
#pragma omp parallel for
for (int step = 0; step < lim / steplenght; step++)
{
for (int i = 0; i < steplenght / 2; i++)
{
int Index0 = steplenght * step + i;
int Index1 = steplenght * step + i + steplenght / 2;
int idx = lim / steplenght * i;
FFT::cpx temp = tempA[Index1] * WN[idx];
tempA[Index1] = tempA[Index0] - temp;
tempA[Index0] = tempA[Index0] + temp;
}
}
}
dst.resize(lim);
for (int i = 0; i < lim; i++)
{
if (opt == -1) dst[i] = tempA[i] / double(lim);
else dst[i] = tempA[i];
}
}
void FFT::fft2d(cpx2d& src, cpx2d& dst, bool inverse, double freq)
{
if (!__is_exp2((int)src.size()) || !__is_exp2((int)src[0].size())) {
FFT::expand(src, FFT::cpx(0, 0));
}
int dim1 = (int)src.size(), dim2 = (int)src[0].size();
FFT::init(dst, dim1, dim2, FFT::cpx(0, 0));
// 沿第一维
#pragma omp parallel for
for (int i = 0; i < dim1; ++i) {
cpx1d temp(dim2);
for (int j = 0; j < dim2; ++j) {
temp[j] = src[i][j];
}
cpx1d temp2;
FFT::fft(temp, temp2, inverse);
for (int j = 0; j < dim2; ++j) {
dst[i][j] = temp2[j];
}
}
// 沿第二维
#pragma omp parallel for
for (int i = 0; i < dim2; ++i) {
cpx1d temp(dim1);
for (int j = 0; j < dim1; ++j) {
temp[j] = dst[j][i];
}
cpx1d temp2;
FFT::fft(temp, temp2, inverse);
for (int j = 0; j < dim1; ++j) {
dst[j][i] = temp2[j];
}
}
}
void FFT::fft3d(cpx3d& src, cpx3d& dst, bool inverse, double freq)
{
if (!__is_exp2((int)src.size()) || !__is_exp2((int)src[0].size()) || !__is_exp2((int)src[0][0].size())) {
FFT::expand(src, FFT::cpx(0, 0));
}
int dim1 = (int)src.size(), dim2 = (int)src[0].size(), dim3 = (int)src[0][0].size();
FFT::init(dst, dim1, dim2, dim3, FFT::cpx(0, 0));
// 沿第一维
#pragma omp parallel for
for (int i = 0; i < dim1; ++i) {
cpx2d temp;
FFT::init(temp, dim2, dim3, FFT::cpx(0, 0));
for (int j = 0; j < dim2; ++j) for (int k = 0; k < dim3; ++k) {
temp[j][k] = src[i][j][k];
}
cpx2d temp2;
FFT::fft2d(temp, temp2, inverse);
for (int j = 0; j < dim2; ++j) for (int k = 0; k < dim3; ++k) {
dst[i][j][k] = temp2[j][k];
}
}
// 沿第二维
#pragma omp parallel for
for (int i = 0; i < dim2; ++i) {
cpx2d temp;
FFT::init(temp, dim1, dim3, FFT::cpx(0, 0));
for (int j = 0; j < dim1; ++j) for (int k = 0; k < dim3; ++k) {
temp[j][k] = dst[j][i][k];
}
cpx2d temp2;
FFT::fft2d(temp, temp2, inverse);
for (int j = 0; j < dim1; ++j) for (int k = 0; k < dim3; ++k) {
dst[j][i][k] = temp2[j][k];
}
}
// 沿第三维
#pragma omp parallel for
for (int i = 0; i < dim3; ++i) {
cpx2d temp;
FFT::init(temp, dim1, dim2, FFT::cpx(0, 0));
for (int j = 0; j < dim1; ++j) for (int k = 0; k < dim2; ++k) {
temp[j][k] = dst[j][k][i];
}
cpx2d temp2;
FFT::fft2d(temp, temp2, inverse);
for (int j = 0; j < dim1; ++j) for (int k = 0; k < dim2; ++k) {
dst[j][k][i] = temp2[j][k];
}
}
}
void FFT::spectrum(cpx2d& src, vector2d& dist)
{
//自动缩放到0-255范围内,并变换象限，将低频移至中间
int dim1 = (int)src.size();
int dim2 = (int)src[0].size();
auto todouble = [](FFT::cpx& cp) {
double norm = cp.real() * cp.real() + cp.imag() * cp.imag();
return log(sqrt(norm) + 1);
};
double max = todouble(src[0][0]);
double min = max;
for (int i = 0; i < dim1; i++) for (int j = 0; j < dim2; ++j)
{
double val = todouble(src[i][j]);;
if (val > max) max = val;
if (val < min) min = val;
}
double scale = 255.0 / (max - min);
FFT::init<double>(dist, dim1, dim2, 0);
for (int i = 0; i < dim1; i++) for (int j = 0; j < dim2; ++j)
{
int target_i = (i + dim1 / 2) % dim1;
int target_j = (j + dim2 / 2) % dim2;
dist[target_i][target_j] = (todouble(src[i][j]) - min) * scale;
}
}
template<typename T>
void FFT::expand(vector<T>& vec1d, T val)
{
vec1d.resize(__expand_exp2((int)vec1d.size()), val);
}
template<typename T>
void FFT::expand(vector<vector<T>>& vec2d, T val)
{
size_t dim1 = vec2d.size(), dim2 = vec2d[0].size();
int new_dim1 = __expand_exp2((int)dim1), new_dim2 = __expand_exp2((int)dim2);
vec2d.resize(new_dim1);
#pragma omp parallel for
for (int i = 0; i < new_dim1; ++i) {
vec2d[i].resize(new_dim2, val);
}
}
template<typename T>
void FFT::expand(vector<vector<vector<T>>>& vec3d, T val)
{
size_t dim1 = vec3d.size(), dim2 = vec3d[0].size(), dim3 = vec3d[0][0].size();
int new_dim1 = __expand_exp2((int)dim1), new_dim2 = __expand_exp2((int)dim2), new_dim3 = __expand_exp2((int)dim3);
vec3d.resize(new_dim1);
#pragma omp parallel for
for (int i = 0; i < new_dim1; ++i) {
vec3d[i].resize(new_dim2);
for (int j = 0; j < new_dim2; ++j) {
vec3d[i][j].resize(new_dim3, val);
}
}
}
template<typename T>
void FFT::dosome(vector<T>& vec1d, function<void(T&)> func)
{
size_t dim1 = vec1d.size();
#pragma omp parallel for
for (int i = 0; i < dim1; ++i) {
func(vec1d[i]);
}
}
template<typename T>
void FFT::dosome(vector<vector<T>>& vec2d, function<void(T&)> func)
{
size_t dim1 = vec2d.size(), dim2 = vec2d[0].size();
#pragma omp parallel for
for (int i = 0; i < dim1; ++i) {
for (int j = 0; j < dim2; ++j) {
func(vec2d[i][j]);
}
}
}
template<typename T>
void FFT::dosome(vector<vector<vector<T>>>& vec3d, function<void(T&)> func)
{
size_t dim1 = vec3d.size(), dim2 = vec3d[0].size(), dim3 = vec3d[0][0].size();
#pragma omp parallel for
for (int i = 0; i < dim1; ++i) {
for (int j = 0; j < dim2; ++j) {
for (int k = 0; k < dim3; ++k) {
func(vec3d[i][j][k]);
}
}
}
}
此处可能存在不合适展示的内容，页面不予展示。您可通过相关编辑功能自查并修改。
如您确认内容无涉及 不当用语 / 纯广告导流 / 暴力 / 低俗色情 / 侵权 / 盗版 / 虚假 / 无价值内容或违法国家有关法律法规的内容，可点击提交进行申诉，我们将尽快为您处理。
马建仓 AI 助手