计算方法实验8：快速傅里叶变换与快速傅里叶逆变换实现对给定函数的 Fourier 分析及重建

千里澄江

已于 2024-06-06 14:38:15 修改

阅读量1.1k

点赞数 21

分类专栏： C/C++ 计算方法实验文章标签：算法 c++ python 图像处理傅里叶分析

于 2024-05-20 14:49:49 首次发布

本文链接：https://blog.csdn.net/weixin_73004416/article/details/139062124

版权

C/C++ 同时被 2 个专栏收录

30 篇文章 1 订阅

订阅专栏

计算方法实验

10 篇文章 0 订阅

订阅专栏

任务

通过快速傅里叶变换与快速傅里叶逆变换实现对给定函数的 Fourier 分析， $t\in[0,1)$ ，函数 $f$ 以及划分数 $n$ 如下：

$f_1( t) = 0. 7\sin ( 2\pi \times 2t) + \sin ( 2\pi \times 5t) , n= 2^4, 2^7；$
$f_2( t) = 0. 7\sin ( 2\pi \times 2t) + \sin ( 2\pi \times 5t) + 0. 3\times random( t)$ ,其中 random(t)为[0,1)区间内的随机数，n= $2^7$ .

对于每个 $f$ 与 $n$ , 程序应输出所得到的快速傅里叶变换后的向量 $g$ , 将 $g$ 的每个分量的实部与虚部分开输出；
对于每个 $f$ , 将取不同的 $n$ 获得的 $g$ 的每个分量模长 $g_i|$ 绘制成图像(为了使图像清晰，如果将不同的 $n$ 绘制到一张图上请绘制折线图),图像横轴为频率(获取方式见附录), 纵轴为 $g_i|$ ；
对于 $f_1$ , 绘制 $f_1$ 离散化后的图像 $f_{1,k})$ , 以及快速傅里叶变换与快速傅里叶逆变换后的最终结果图像。不同的 $n$ 需绘制在不同的图像上，绘制方式为折线图；
对于 $f_2$ ,绘制 $f_2$ 离散化后的图像 $f_{2,k})$ , 快速傅里叶变换与快速傅里叶逆变换后的最终结果图像，以及快速傅里叶变换后取频率域前25%的系数进行快速傅里叶逆变换所得最终结果图像，绘制方式为折线图。

算法

1.FFT

$n\leftarrow length[f]$

$\textbf{if }n= = 1$

$\quad\quad\textbf{then return }f$

$\textbf{end if}$

$\omega_n\leftarrow e^{\boldsymbol{i}2\pi/n}$

$\omega\leftarrow1$

$\mathbf{f}^{0}\leftarrow(f_{0},f_{2},\ldots,f_{n-2})$

$\mathbf{f}^{1}\leftarrow(f_{1},f_{3},\ldots,f_{n-1})$

$\mathbf{g}^0\leftarrow$ FFT $(\mathbf{f}^0)$

$\mathbf{g}^1\leftarrow$ FFT $(\mathbf{f}^1)$

$\textbf{for }k\leftarrow 0$ to $1\textbf{do}$

$\quad\quad\mathbf{g_k}\leftarrow\mathbf{g}_k^0+\omega\mathbf{g}_k^1$

$\quad\quad\mathbf{g_{k+n/2}}\leftarrow\mathbf{g}_k^0-\omega\mathbf{g}_k^1$

$\quad\quad\omega\leftarrow\omega\omega_n$

$\textbf{end for}$

$\textbf{return g}$

2.IFFT

方法1

$n\leftarrow length[f]$

$\textbf{if }n= = 1$

$\quad\quad\textbf{then return }f$

$\textbf{end if}$

$\omega_n\leftarrow e^{\boldsymbol{-i}2\pi/n}$

$\omega\leftarrow1$

$\mathbf{f}^{0}\leftarrow(f_{0},f_{2},\ldots,f_{n-2})$

$\mathbf{f}^{1}\leftarrow(f_{1},f_{3},\ldots,f_{n-1})$

$\mathbf{g}^0\leftarrow$ IFFT $(\mathbf{f}^0)$

$\mathbf{g}^1\leftarrow$ IFFT $(\mathbf{f}^1)$

$\textbf{for }k\leftarrow 0$ to $1\textbf{do}$

$\quad\quad\mathbf{g_k}\leftarrow\mathbf{g}_k^0+\omega\mathbf{g}_k^1$

$\quad\quad\mathbf{g_{k+n/2}}\leftarrow\mathbf{g}_k^0-\omega\mathbf{g}_k^1$

$\quad\quad\omega\leftarrow\omega\omega_n$

$\textbf{end for}$

$\mathbf{g} = \mathbf{g}/2$

$\textbf{return g}$

方法2

$\mathbf{f_{temp}}\leftarrow(f_{0},f_{n-1},f_{n-2}\ldots,f_{2})$
$\mathbf{g}\leftarrow$ FFT $(\mathbf{f_{temp}})$
$\mathbf{g} = \mathbf{g}/n$
$\textbf{return g}$

代码

C++实现FFT与IFFT

#include<bits/stdc++.h>

using namespace std;

const double PI = acos(-1);
random_device rd;
mt19937 gen(rd());
uniform_real_distribution<> dis(0, 1);

complex<double> f1(double t);
complex<double> f2(double t);
void fft(vector<complex<double>> &a);
void ifft(vector<complex<double>> &a);

int main() {
    ofstream fft_file;
    ofstream discrete;
    ofstream ifft_file;
    cout << "Function 1:\n";
    for (int n : {16, 128}) {
        discrete.open(n == 16 ? "discrete_f1(1).txt" : "discrete_f1(2).txt");
        vector<complex<double>> g1(n);
        for (int i = 0; i < n; i++) {
            g1[i] = f1(i / (double)n);
            discrete << i / (double)n << ' ' << g1[i].real() << endl;
        }
        discrete.close();

        fft(g1);
        
        fft_file.open(n == 16 ? "fft_f1(1).txt" : "fft_f1(2).txt");
        for (int i = 0; i < n; i++) {
            if (i < n / 2)
                fft_file << i + 1 << ' ' << abs(g1[i]) << endl;
            else
                fft_file << i - n << ' ' << abs(g1[i]) << endl;
        }
            
        fft_file.close();
        cout << "FFT of function for n = " << n << endl;
        cout << "Real part:\n";
        for (const auto& x : g1)
            cout << x.real() << ' ';
        cout << endl;
        cout << "Imaginary part:\n";
        for (const auto& x : g1)
            cout << x.imag() << ' ';
        cout << endl << endl;

        vector<complex<double>> inv1 = g1;
        ifft(inv1);
        ifft_file.open(n == 16 ? "ifft_f1(1).txt" : "ifft_f1(2).txt");
        for (int i = 0; i < n; i++)
            ifft_file << i / (double)n << ' ' << inv1[i].real() << endl;
        ifft_file.close();
    }
    fft_file.close();


    cout << "Function 2:\n";
    vector<complex<double>> g2(128);
    discrete.open("discrete_f2.txt");
    for (int i = 0; i < 128; i++) {
        g2[i] = f2((double)i / 128);
        discrete << i / 128.0 << ' ' << g2[i].real() << endl;
    }
    discrete.close();

    fft(g2);

    fft_file.open("fft_f2.txt");
    for (int i = 0; i < 128; i++){
        if (i < 64)
            fft_file << i + 1 << ' ' << abs(g2[i]) << endl;
        else
            fft_file << i - 128 << ' ' << abs(g2[i]) << endl;
    }
    cout << "FFT of function for n = " << 128 << endl;
    cout << "Real part:\n";
    for (const auto& x : g2)
        cout << x.real() << ' ';
    cout << endl;
    cout << "Imaginary part:\n";
    for (const auto& x : g2)
        cout << x.imag() << ' ';
    cout << endl;
    
    vector<complex<double>> inv2 = g2;
    ifft(inv2);
    ifft_file.open("ifft_f2(1).txt");
    for (int i = 0; i < 128; i++)
        ifft_file << (double)i / 128 << ' ' << inv2[i].real() << endl;
    ifft_file.close();

    vector<complex<double>> inv3(g2.begin(), g2.begin() + 32);
    inv3.resize(128, complex<double>(0, 0));
    ifft(inv3);
    ifft_file.open("ifft_f2(2).txt");
    for (int i = 0; i < 128; i++)
        ifft_file << (double)i / 128 << ' ' << inv3[i].real() << endl;
    ifft_file.close();
}

complex<double> f1(double t) {
    return 0.7 * sin(2 * PI * 2 * t) + sin(2 * PI * 5 * t);
}

complex<double> f2(double t) {
    return 0.7 * sin(2 * PI * 2 * t) + sin(2 * PI * 5 * t) + 0.3 * dis(gen);
}

void fft(vector<complex<double>> &a) {
    int n = a.size();
    if (n == 1) 
        return;

    vector<complex<double>> a0(n / 2), a1(n / 2);
    for (int i = 0; 2 * i < n; i++) {
        a0[i] = a[2 * i];
        a1[i] = a[2 * i + 1];
    }

    fft(a0);
    fft(a1);
    
    double ang = 2 * PI / n;
    complex<double> w(1), wn(cos(ang), sin(ang));
    for (int i = 0; 2 * i < n; i++) {
        a[i] = a0[i] + w * a1[i];
        a[i + n / 2] = a0[i] - w * a1[i];
        w *= wn;
    }
}

/*void ifft(vector<complex<double>> &a) {
    int n = a.size();
    reverse(a.begin() + 1, a.end());

    fft(a);

    for (auto& x : a)
        x /= n;
}*/

void ifft(vector<complex<double>> &a) {
    int n = a.size();
    if (n == 1)
        return;

    vector<complex<double>> a0(n / 2), a1(n / 2);
    for (int i = 0; 2 * i < n; i++) {
        a0[i] = a[2 * i];
        a1[i] = a[2 * i + 1];
    }

    ifft(a0);
    ifft(a1);

    double ang = -2 * PI / n; // negative sign
    complex<double> w(1), wn(cos(ang), sin(ang));
    for (int i = 0; 2 * i < n; i++) {
        a[i] = a0[i] + w * a1[i];
        a[i + n / 2] = a0[i] - w * a1[i];
        w *= wn;
    }

    for (auto& x : a)
        x /= 2;
}

python实现结果可视化

f1 FFT

import matplotlib.pyplot as plt

# 读取数据
def read_data(file_path):
    with open(file_path, 'r') as file:
        data = file.readlines()
    x = []
    y = []
    for line in data:
        values = line.split()
        x.append(float(values[0]))
        y.append(float(values[1]))
    return x, y

# 从文件中读取数据
x1, y1 = read_data('fft_f1(1).txt')
x2, y2 = read_data('fft_f1(2).txt')

# 绘制图像
plt.figure(figsize=(12, 6))

plt.plot(x1, y1, marker='o', color='blue', label='n=16')
plt.plot(x2, y2, marker='o', color='red', label='n=128')

plt.title('FFT of f1')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()

plt.tight_layout()
plt.show()

f2 FFT

import matplotlib.pyplot as plt

def read_data(file_path):
    with open(file_path, 'r') as file:
        data = file.readlines()
    x = []
    y = []
    for line in data:
        values = line.split()
        x.append(float(values[0]))
        y.append(float(values[1]))
    return x, y

x, y = read_data('fft_f2.txt')

# 绘制图像
plt.plot(x, y, marker='o', label='FFT of f2')

plt.title('FFT of f2')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()

plt.show()

Discrete and IFFT of f1(1)

import matplotlib.pyplot as plt

# 读取数据
def read_data(file_path):
    with open(file_path, 'r') as file:
        data = file.readlines()
    x = []
    y = []
    for line in data:
        values = line.split()
        x.append(float(values[0]))
        y.append(float(values[1]))
    return x, y

# 从文件中读取数据
x_discrete, y_discrete = read_data('discrete_f1(1).txt')
x_ifft, y_ifft = read_data('ifft_f1(1).txt')

# 绘制图像
plt.figure(figsize=(12, 6))

plt.plot(x_discrete, y_discrete, marker='o', color='blue', label='Discrete f1')
plt.plot(x_ifft, y_ifft, marker='o', color='red', label='IFFT of f1')

plt.title('Discrete and IFFT of f1')
plt.xlabel('x')
plt.ylabel('f1')
plt.legend()

plt.tight_layout()
plt.show()

Discrete and IFFT of f1(2)

import matplotlib.pyplot as plt

# 读取数据
def read_data(file_path):
    with open(file_path, 'r') as file:
        data = file.readlines()
    x = []
    y = []
    for line in data:
        values = line.split()
        x.append(float(values[0]))
        y.append(float(values[1]))
    return x, y

# 从文件中读取数据
x_discrete, y_discrete = read_data('discrete_f1(2).txt')
x_ifft, y_ifft = read_data('ifft_f1(2).txt')

# 绘制图像
plt.figure(figsize=(12, 6))

plt.plot(x_discrete, y_discrete, marker='o', color='blue', label='Discrete f1')
plt.plot(x_ifft, y_ifft, marker='o', color='red', label='IFFT of f1')

plt.title('Discrete and IFFT of f1')
plt.xlabel('x')
plt.ylabel('f1')
plt.legend()

plt.tight_layout()
plt.show()

Discrete and IFFT of f2(1)

import matplotlib.pyplot as plt

# 读取数据
def read_data(file_path):
    with open(file_path, 'r') as file:
        data = file.readlines()
    x = []
    y = []
    for line in data:
        values = line.split()
        x.append(float(values[0]))
        y.append(float(values[1]))
    return x, y

# 从文件中读取数据
x_discrete, y_discrete = read_data('discrete_f2.txt')
x_ifft, y_ifft = read_data('ifft_f2(1).txt')

# 绘制图像
plt.figure(figsize=(12, 6))

plt.plot(x_discrete, y_discrete, marker='o', color='blue', label='Discrete f2')
plt.plot(x_ifft, y_ifft, marker='o', color='red', label='IFFT of f2')

plt.title('Discrete and IFFT of f2')
plt.xlabel('x')
plt.ylabel('f1')
plt.legend()

plt.tight_layout()
plt.show()

Discrete and IFFT of f2(2)

import matplotlib.pyplot as plt

# 读取数据
def read_data(file_path):
    with open(file_path, 'r') as file:
        data = file.readlines()
    x = []
    y = []
    for line in data:
        values = line.split()
        x.append(float(values[0]))
        y.append(float(values[1]))
    return x, y

# 从文件中读取数据
x_discrete, y_discrete = read_data('discrete_f2.txt')
x_ifft, y_ifft = read_data('ifft_f2(2).txt')

# 绘制图像
plt.figure(figsize=(12, 6))

plt.plot(x_discrete, y_discrete, marker='o', color='blue', label='Discrete f2')
plt.plot(x_ifft, y_ifft, marker='o', color='red', label='IFFT of f2')

plt.title('Discrete and IFFT of f2')
plt.xlabel('x')
plt.ylabel('f2')
plt.legend()

plt.tight_layout()
plt.show()

python直接实现FFT和IFFT及结果的可视化

import numpy as np
import matplotlib.pyplot as plt

# 定义函数
def f1(t):
    return 0.7 * np.sin(2 * np.pi * 2 * t) + np.sin(2 * np.pi * 5 * t)

# 设置采样率和时间范围
sample_rates = [128]  # 只使用抽样频率为128的情况
time_duration = 1  # 总时间范围

# 对每个采样率进行FFT并绘图
for sample_rate in sample_rates:
    # 生成采样点
    time_range = np.arange(0, time_duration, 1/sample_rate)
    samples = f1(time_range)

    # 计算FFT
    fft_result = np.fft.fft(samples)
    fft_result = np.fft.fftshift(fft_result)  # Shift zero freq to center

    # 计算频率
    freqs = np.fft.fftfreq(len(samples), d=1/sample_rate)
    freqs = np.fft.fftshift(freqs)  # Shift zero freq to center

    # 计算IFFT
    ifft_result = np.fft.ifft(np.fft.ifftshift(fft_result))

    # 创建一个新的图形窗口
    plt.figure(figsize=(10, 5))

    # 在同一图上绘制f1的离散点折线图和IFFT图像
    plt.plot(time_range, samples, label='f1 discrete points')
    plt.plot(time_range, ifft_result.real, label='IFFT result')
    plt.title(f'f1 discrete points and IFFT result with sample rate {sample_rate}')
    plt.xlabel('Time (s)')
    plt.ylabel('Amplitude')
    plt.grid(True) # 显示网格
    plt.legend()

    # 创建一个新的图形窗口
    plt.figure(figsize=(10, 5))

    # 绘制FFT图像
    plt.plot(freqs, np.abs(fft_result), label='FFT result')
    plt.title(f'FFT of f1 with sample rate {sample_rate}')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Amplitude')
    plt.grid(True)
    plt.legend()

    # 显示图形
    plt.show()