北邮22信通一枚~
跟随课程进度更新北邮信通院DSP的笔记、代码和文章,欢迎关注~
获取更多文章,请访问专栏:
目录
一、实验内容
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算法差异
运行效果