北邮22级信通院DSP：实验三（2）：C++、MATLAB对DFT和FFT的实现、DFT和FFT性能比较、DFDM拓展实验N=64和N=128点算法差异（MATLAB绘图）

最新推荐文章于 2024-05-30 20:15:33 发布

青山入墨雨如画

最新推荐文章于 2024-05-30 20:15:33 发布

阅读量973

点赞数 21

分类专栏：北邮22级信通院DSP 文章标签：算法 matlab dsp开发 FFT DFT c++ 开发语言

本文链接：https://blog.csdn.net/bc202205/article/details/138861307

版权

北邮22级信通院DSP 专栏收录该内容

6 篇文章 1 订阅

订阅专栏

北邮22信通一枚~

跟随课程进度更新北邮信通院DSP的笔记、代码和文章，欢迎关注~

获取更多文章，请访问专栏：

北邮22级信通院DSP_青山入墨雨如画的博客-CSDN博客

一、实验内容

1.1实验准备

实验需要使用C++、C环境在OFDM场景下实现FFT与IFFT算法，然后利用MATLAB对频谱进行绘制和分析。

1.2实验步骤

（1）根据FFT、IFFT定义以及DFT定义分别编写对应变换的函数；

（2）生成2^n个子载波（其中n>6）,利用MATLAB画出频域示意图，标注频域间隔；

（3）对生成的2^n个子载波进行IFFT变换，输出变换结果，并在MATLAB中绘图展示；

（4）将调制后的信号s（t）进行FFT变换，输出变换结果，并在MATLAB中绘图展示；

（5）上述步骤分别计算使用FFT以及DFT的运行结果，比较两种算法的性能；

（6）利用MATLAB自带的FFT函数进行验证，检查算法结果的差异。

1.3拓展实验

对实验结果进行自主分析，比较N=64和N=128算法差异。

二、实验要求

建议使用软件：MATLAB，visual studio（能编译运行C、C++软件均可）

三、C++算法及代码

3.1DFT算法

void dft(vector<Complex>&a)
{
    complex<double>W, Wk;
    Wk = 1;
    int n = a.size();
    W.real(cos(PI * 2 / n));
    W.imag(-sin(PI * 2 / n));
    complex<double>sum = 0;
    for (int i = 0; i < n; i++)
    {
        sum = 0;
        for (int j = 0; j < n; j++)
        {
            sum += a[j] * pow(Wk, j);
        }
        Wk *= W;
        a[i] = sum;
    }
    return;
}

3.2FFT算法

void fft(vector<Complex>& a, bool inverse = false)
{
    int n = a.size();                         // 获取输入向量的大小
    if (n <= 1) return;                       // 如果输入向量大小为1或0，则直接返回，无需变换

    // 分别定义偶数项和奇数项的向量
    vector<Complex> even(n / 2), odd(n / 2);
    for (int i = 0; i < n / 2; ++i)
    {
        even[i] = a[i * 2];                   // 偶数项
        odd[i] = a[i * 2 + 1];                // 奇数项
    }

    fft(even, inverse);                       // 对偶数项进行递归FFT变换
    fft(odd, inverse);                        // 对奇数项进行递归FFT变换

    // 计算旋转因子的角度
    double angle = (inverse ? 2 : -2) * PI / n;       //false的话是-2，正变换;true的话是2，逆变换
    Complex w(1), wn(cos(angle), sin(angle));         // 初始化旋转因子
    for (int i = 0; i < n / 2; ++i)
    {
        a[i] = even[i] + w * odd[i];                  // 计算上半部分的结果
        a[i + n / 2] = even[i] - w * odd[i];          // 计算下半部分的结果
        if (inverse)
        {                                             // 如果是逆变换，则需要除以2
            a[i] /= 2;
            a[i + n / 2] /= 2;
        }
        w *= wn; // 更新旋转因子
    }
}

3.3FFT算法讲解

北邮22级信通院DSP：实验三（1）：FFT变换、IFFT变换（附每步8点变换蝶形图）保姆级讲解+用C++程序实现复数域的FFT变换和IFFT变换+C++中的chrono头文件讲解-CSDN博客

3.4完整代码

#include<iostream>                                        // 输入输出流头文件
#include <complex>                                        // 包含复数类型的头文件
#include <vector>                                         // 包含向量容器的头文件
#include <cmath>                                          // 包含数学函数的头文件
#include <chrono>                                         // 用于处理程序运行时间的头文件
#include <functional>                                     // 引入 <functional> 头文件用于使用 std::function方法
using namespace std;

//用反三角函数定义常数PI
const double PI = acos(-1.0);

// 定义复数类型
typedef complex<double> Complex;

// FFT算法函数，参数 a 是输入的复数向量，
// inverse 表示是否进行逆变换，默认为 false
void fft(vector<Complex>& a, bool inverse = false)
{
    int n = a.size();                         // 获取输入向量的大小
    if (n <= 1) return;                       // 如果输入向量大小为1或0，则直接返回，无需变换

    // 分别定义偶数项和奇数项的向量
    vector<Complex> even(n / 2), odd(n / 2);
    for (int i = 0; i < n / 2; ++i)
    {
        even[i] = a[i * 2];                   // 偶数项
        odd[i] = a[i * 2 + 1];                // 奇数项
    }

    fft(even, inverse);                       // 对偶数项进行递归FFT变换
    fft(odd, inverse);                        // 对奇数项进行递归FFT变换

    // 计算旋转因子的角度
    double angle = (inverse ? 2 : -2) * PI / n;       //false的话是-2，正变换;true的话是2，逆变换
    Complex w(1), wn(cos(angle), sin(angle));         // 初始化旋转因子
    for (int i = 0; i < n / 2; ++i)
    {
        a[i] = even[i] + w * odd[i];                  // 计算上半部分的结果
        a[i + n / 2] = even[i] - w * odd[i];          // 计算下半部分的结果
        if (inverse)
        {                                             // 如果是逆变换，则需要除以2
            a[i] /= 2;
            a[i + n / 2] /= 2;
        }
        w *= wn; // 更新旋转因子
    }
}

// IFFT算法函数，参数 a 是输入的复数向量
void ifft(vector<Complex>& a)
{
    fft(a, true); // 调用 FFT 函数，inverse 参数设置为 true，进行逆变换
}


void dft(vector<Complex>&a)
{
    complex<double>W, Wk;
    Wk = 1;
    int n = a.size();
    W.real(cos(PI * 2 / n));
    W.imag(-sin(PI * 2 / n));
    complex<double>sum = 0;
    for (int i = 0; i < n; i++)
    {
        sum = 0;
        for (int j = 0; j < n; j++)
        {
            sum += a[j] * pow(Wk, j);
        }
        Wk *= W;
        a[i] = sum;
    }
    return;
}
// 计算运行时间
template<typename Func>
void measureTime(Func func, const string& msg) {
    auto start = chrono::high_resolution_clock::now();
    func(); // 执行传入的函数
    auto end = chrono::high_resolution_clock::now();
    chrono::duration<double> duration = end - start;
    cout << msg << "运行时间：" << duration.count() * 1000 << " ms" << endl;
}

// 计算内存占用
template<typename T>
size_t memoryUsage(const vector<T>& vec) {
    return vec.size() * sizeof(T);
}

int main()
{
    system("color 0A");
    // 输入序列
    vector<Complex> input = { (1,2),(2,3),(3,4),(4,5),(5,6),(6,7),(7,8),(8,9),(9,1) };
    // 执行FFT变换
    fft(input);
    // 打印FFT变换结果
    cout << "FFT 变换结果：" << endl << endl;
    for (const auto& val : input)
    {
        cout << val << endl;
    }
    cout << endl;

    // 执行IFFT变换
    ifft(input);

    // 打印IFFT变换结果
    cout << "IFFT 变换结果：" << endl << endl;
    for (const auto& val : input)
    {
        cout << val << endl;
    }
    cout << endl;

    // 计算 IFFT 变换的运行时间
    measureTime([&input]() {ifft(input); }, "FFT");
    measureTime([&input]() {dft(input); }, "DFT");
    cout << endl;
    // 计算内存占用
    cout << "FFT内存占用：" << memoryUsage(input) << " bit" << endl;
    cout << "DFT内存占用：" << memoryUsage(input) << " bit" << endl;
    return 0;
}

四、MATLAB算法及代码

4.1第一种写法

代码部分

clear all

%%
%FFT与DFT对比实验
n = 10;
N = 2^n; % 变换点数 N
t = 0:1:N-1; % 采样点 t
x = sin(3.*t); % 离散信号 x
% 测试运行时间
tic
use_mydft = mydft(x); % For 循环的 DFT 变换
toc
%历时 0.121736 秒。

tic
use_myfft = myfft(x, N); % 基 2 时间抽选法的 FFT 变换
toc
%历时 0.015503 秒。

tic
use_matlabfft = fft(x); % MATLAB 的 FFT 变换
toc
%历时 0.003128 秒。

%%
%OFDM拓展实验
%%

% 1、OFDM 子载波 IFFT 变换并绘制波形图
%{
首先生成 2^n 个 OFDM 子载波，n 取 10，
为简化实验，仅考虑前 6 个子载波进行变换。
使用 IFFT 变换公式求得载波变换后的结果。
%}
Fs = 1000; % 采样频率
N = 2^10; % 子载波数
T = N/Fs; % 信号绘制为一个周期的长度
t = 0:1/Fs:T-1/Fs; 
num = 6; % 仅考虑 6 个子载波的变换
phase = repmat((1+1i)/sqrt(2),1,T*Fs); % 初始相位
y(num,T*Fs-1) = 0; % 先分配空间加快程序运行速度
for k=0:num-1
 for n=0:T*Fs-1
 y(k+1,n+1)=phase(n+1)*exp(1i*2*pi*k*n/N); % IFFT 变换
 end
end

%其次，绘制 IFFT 信号的实部和虚部。
figure(1)
% 绘制 IFFT 信号的实部
subplot(2,1,1)
for k=0:num-1
 plot(0:T*Fs-1, real(y(k+1,:)))
 hold on;
end
hold off;
title('IFFT 后信号波形图（实部）') % 标题
xlabel('时间') % 横轴名称
ylabel('载波') % 纵轴名称
set(gca,'XTick',0:0:0); % 去掉横轴刻度
set(gca,'YTick',-1:1:1); % 去掉纵轴刻度
axis([0 N -1 1]) % 横轴纵轴坐标范围
% 绘制 IFFT 信号的虚部
subplot(2,1,2)
for k=0:num-1
 plot(0:T*Fs-1, imag(y(k+1,:)))
 hold on;
end
hold off;
title('IFFT 后信号波形图（虚部）') % 标题
xlabel('时间') % 横轴名称
ylabel('载波') % 纵轴名称
set(gca,'XTick',0:0:0); % 去掉横轴刻度
set(gca,'YTick',-1:1:1); % 去掉纵轴刻度
axis([0 N -1 1]) % 横轴纵轴坐标范围

%%
%2、调制后信号 FFT 变换并绘制频谱图
%首先，对上述信号求 FFT 变换
zero_num = 20; % 补零位数
y_zero = [y,zeros(num,zero_num*1024)]; % 对 y 末尾补 0
f = (0:(zero_num+1)*T*Fs-1)/T/(zero_num+1)-Fs/2;
for k=0:num-1
 Y(k+1,:)=abs(fftshift(fft(y_zero(k+1,:))))/N; % FFT 变换
end

%其次，绘制信号频谱图
figure(2)
for k=0:num-1
 plot(f,Y(k+1,:))
 hold on;
end
hold off;
grid on; % 网格
title('OFDM 信号频谱图') % 标题
xlim([-10,10]);
xlabel('频率/Hz') % 横轴名称
ylabel('幅度') % 纵轴名称

%%
%附录

%For 循环的 DFT 变换函数 mydft.m 文件
function [X] = mydft(x)
% For 循环写 DFT 变换
for i = 0:length(x)-1
 y = 0;
 for j = 0:length(x)-1
 y = y + x(j+1)*exp(-1j*2*pi*i*j/length(x));
 end
 X(i+1) = y;
end
end

%基 2 时间抽选法的 FFT 变换函数 myfft.m 文件
function X = myfft(x,M)
% 基 2 时间抽选法写 FFT 变换，M 必须为 2 的 n 次方，可以不输入 M
x_len = length(x);
% 根据是否输入 M 将 x 补到 2 的 n 次方长度
if nargin<2
 % 如果参数少于 2
 n = 1;
while(x_len > 2^n)
 n = n+1;
 end
 x_zero = zeros(1,2^n);
else
 % 如果参数大于等于 2
 x_zero = zeros(1,M);
 n = log2(M);
end
x_zero(1:x_len) = x;
x_zero_len = length(x_zero);
% 对位置进行二进制取反后转为 10 进制
x = bin2dec(fliplr(dec2bin(0:x_zero_len-1,n)))+1;
x1 = x_zero(x);
% 碟式变换的级数
for m1 = 1:n
 % 基 2 时间抽选法，数值的间隔
 space = 2^(m1-1);
 for m2 = 0:space-1
 % 奇数项的旋转因子的指数
 k = 2^(n-m1)*m2;
 for m3 = m2+1:2^m1:x_zero_len
 % 奇数项旋转因子
 W = exp(-1j*2*pi*k/x_zero_len);
 % 先存一下变换后的 x1(m3)
 y = x1(m3)+x1(m3+space)*W;
 % 重复对 x1 调用即可，无需新的变量
 x1(m3+space) = x1(m3)-x1(m3+space)*W;
 x1(m3) = y;
 end
 end
end
X = x1;
end

运行效果

4.2第二种写法

代码部分

% 步骤一：生成2^n个子载波并绘制频域示意图

% 参数设置
n = 6; % 子载波数量的指数，这里设置为6，表示生成64个子载波
Fs = 1000; % 采样频率
f_spacing = Fs / 64; % 子载波频率间隔

% 生成子载波
t = 0:1/Fs:1-1/Fs; % 时间向量
subcarriers = zeros(length(t), 2^n); % 存储子载波信号
for i = 1:2^n
    subcarriers(:,i) = cos(2*pi*f_spacing*(i-1)*t); % 生成单个子载波信号
end

% 绘制频域示意图
figure;
for i = 1:2^n
    subplot(2^n/2, 2, i);
    plot(abs(fft(subcarriers(:,i)))); % 绘制子载波频谱
    title(['子载波 ', num2str(i)]);
    xlabel('频率（Hz）');
    ylabel('幅度');
end
sgtitle('2^n个子载波的频域示意图');

% 步骤二：对生成的2^n个子载波进行IFFT变换并绘图展示

% IFFT变换
time_domain_signal = ifft(subcarriers);

% 绘制IFFT变换结果
figure;
for i = 1:2^n
    subplot(2^n/2, 2, i);
    plot(abs(time_domain_signal(:,i))); % 绘制IFFT变换结果的幅度
    title(['子载波 ', num2str(i)]);
    xlabel('时间');
    ylabel('幅度');
end
sgtitle('2^n个子载波的IFFT变换结果');

% 步骤三：对调制后的信号进行FFT变换并绘图展示

% 生成调制后的信号
modulated_signal = sum(subcarriers, 2);

% FFT变换
modulated_spectrum = fft(modulated_signal);

% 绘制FFT变换结果
figure;
plot(abs(modulated_spectrum)); % 绘制FFT变换结果的幅度
title('调制后的信号的FFT变换结果');
xlabel('频率（Hz）');
ylabel('幅度');

% 步骤四：利用MATLAB自带的FFT函数进行验证

% 使用MATLAB自带的FFT函数进行验证
matlab_fft_result = fft(modulated_signal);

% 检查与自己编写的C/C++算法结果的差异
diff = sum(abs(matlab_fft_result - modulated_spectrum));
disp(['与自己编写的C/C++算法结果的差异为：', num2str(diff)]);

% 步骤五：扩展实验 - 对实验结果进行自主分析

% 根据实验结果进行自主分析，比较N = 64与N = 128算法差异

运行效果

青山入墨雨如画

关注

21
点赞
踩
12

收藏

觉得还不错? 一键收藏
打赏
0
评论
北邮22级信通院DSP：实验三（2）：C++、MATLAB对DFT和FFT的实现、DFT和FFT性能比较、DFDM拓展实验N=64和N=128点算法差异（MATLAB绘图）

（1）根据FFT、IFFT定义以及DFT定义分别编写对应变换的函数；（2）生成2^n个子载波（其中n>6）,利用MATLAB画出频域示意图，标注频域间隔；（3）对生成的2^n个子载波进行IFFT变换，输出变换结果，并在MATLAB中绘图展示；（4）将调制后的信号s（t）进行FFT变换，输出变换结果，并在MATLAB中绘图展示；（5）上述步骤分别计算使用FFT以及DFT的运行结果，比较两种算法的性能；（6）利用MATLAB自带的FFT函数进行验证，检查算法结果的差异。
复制链接

扫一扫