频域分析的核心——从时域采样到频域离散谱
离散傅里叶变换（DFT）将有限长离散信号变换到离散频域，是数字信号处理中最重要的变换。它使得我们可以在频域分析和处理信号。
DFT是DTFT在频率上的等间隔采样：
这意味着DFT的频率分辨率为 Δf = fs/N Hz，增加N可提高频率分辨率。
| 性质 | 时域 | 频域(DFT) |
|---|---|---|
| 线性 | a·x₁[n]+b·x₂[n] | a·X₁[k]+b·X₂[k] |
| 循环移位 | x[(n-m)_N] | W_N^(km)·X[k] |
| 循环卷积 | x[n] ⊛ h[n] | X[k]·H[k] |
| Parseval定理 | Σ|x[n]|² | (1/N)Σ|X[k]|² |
| 对称性 | 实序列x[n] | X[k] = X*[N-k]（共轭对称） |
对无限长信号截断等价于加窗，窗函数的旁瓣会导致频谱泄漏：
| 窗类型 | 主瓣宽度 | 旁瓣衰减(dB) | 应用场景 |
|---|---|---|---|
| 矩形窗 | 4π/N | -13 | 频率分辨率优先 |
| Hanning窗 | 8π/N | -31 | 通用频谱分析 |
| Hamming窗 | 8π/N | -41 | 语音处理 |
| Blackman窗 | 12π/N | -57 | 高动态范围 |
| Kaiser窗 | 可调 | 可调 | 灵活设计 |
实现一个N点DFT的直接计算引擎，使用CORDIC旋转因子生成器：
//=============================================
// dsp_dft_engine.v
// N点DFT直接计算引擎
// X[k] = Σ x[n] * W_N^(nk)
// W_N = e^(-j2π/N) = cos(2πnk/N) - j*sin(2πnk/N)
//=============================================
module dsp_dft_engine #(
parameter N = 16, // DFT点数
parameter DATA_WIDTH = 16, // 输入数据位宽
parameter OUT_WIDTH = 32 // 输出位宽
)(
input wire clk,
input wire rst_n,
input wire start,
input wire signed [DATA_WIDTH-1:0] x_in [0:N-1], // 输入采样
output reg signed [OUT_WIDTH-1:0] X_re [0:N-1], // DFT实部输出
output reg signed [OUT_WIDTH-1:0] X_im [0:N-1], // DFT虚部输出
output reg [OUT_WIDTH-1:0] X_mag [0:N-1], // 幅度谱
output reg done
);
// 旋转因子查找表
reg signed [DATA_WIDTH-1:0] wn_cos_table [0:N*N-1];
reg signed [DATA_WIDTH-1:0] wn_sin_table [0:N*N-1];
initial begin
integer idx, n_idx, k_idx;
for (k_idx = 0; k_idx < N; k_idx = k_idx + 1) begin
for (n_idx = 0; n_idx < N; n_idx = n_idx + 1) begin
idx = k_idx * N + n_idx;
wn_cos_table[idx] = $rtoi(
(2.0**(DATA_WIDTH-1)-1) * $cos(2.0*3.14159265*n_idx*k_idx/N)
);
wn_sin_table[idx] = $rtoi(
(2.0**(DATA_WIDTH-1)-1) * $sin(2.0*3.14159265*n_idx*k_idx/N)
);
end
end
end
// 计算状态机
localparam IDLE = 3'd0;
localparam COMPUTE = 3'd1;
localparam MAGNITUDE = 3'd2;
localparam FINISH = 3'd3;
reg [2:0] state;
reg [$clog2(N)-1:0] k_idx, n_idx;
reg signed [OUT_WIDTH-1:0] re_acc, im_acc;
integer i;
always @(posedge clk or negedge rst_n) begin
if (!rst_n) begin
state <= IDLE;
k_idx <= 0; n_idx <= 0;
done <= 1'b0;
for (i = 0; i < N; i = i + 1) begin
X_re[i] <= 0; X_im[i] <= 0; X_mag[i] <= 0;
end
end else begin
case (state)
IDLE: begin
done <= 1'b0;
if (start) begin
k_idx <= 0;
n_idx <= 0;
re_acc <= 0;
im_acc <= 0;
state <= COMPUTE;
end
end
COMPUTE: begin
// X[k] += x[n] * (cos - j*sin)
re_acc <= re_acc + (x_in[n_idx] * wn_cos_table[k_idx*N+n_idx]);
im_acc <= im_acc - (x_in[n_idx] * wn_sin_table[k_idx*N+n_idx]);
if (n_idx == N - 1) begin
X_re[k_idx] <= re_acc;
X_im[k_idx] <= im_acc;
re_acc <= 0;
im_acc <= 0;
n_idx <= 0;
if (k_idx == N - 1) begin
k_idx <= 0;
state <= MAGNITUDE;
end else begin
k_idx <= k_idx + 1'b1;
end
end else begin
n_idx <= n_idx + 1'b1;
end
end
MAGNITUDE: begin
// 幅度谱：|X[k]| ≈ |X_re| + |X_im|（近似计算）
X_mag[k_idx] <= (X_re[k_idx][OUT_WIDTH-1] ? -X_re[k_idx] : X_re[k_idx])
+ (X_im[k_idx][OUT_WIDTH-1] ? -X_im[k_idx] : X_im[k_idx]);
if (k_idx == N - 1) begin
state <= FINISH;
end else begin
k_idx <= k_idx + 1'b1;
end
end
FINISH: begin
done <= 1'b1;
state <= IDLE;
end
endcase
end
end
endmodule//=============================================
// dsp_window_generator.v
// 可配置窗函数生成器
// 支持：矩形/Hanning/Hamming/Blackman窗
//=============================================
module dsp_window_generator #(
parameter N = 256,
parameter DATA_WIDTH = 16
)(
input wire clk,
input wire rst_n,
input wire [1:0] window_type, // 00=Rect 01=Hanning 10=Hamming 11=Blackman
input wire en,
input wire [$clog2(N)-1:0] index, // 当前采样索引 n
output reg [DATA_WIDTH-1:0] window_out // 窗函数值
);
// 窗函数系数查找表（归一化到DATA_WIDTH位宽）
reg [DATA_WIDTH-1:0] hanning_lut [0:N-1];
reg [DATA_WIDTH-1:0] hamming_lut [0:N-1];
reg [DATA_WIDTH-1:0] blackman_lut [0:N-1];
initial begin
integer i;
for (i = 0; i < N; i = i + 1) begin
// Hanning: 0.5 * (1 - cos(2πn/(N-1)))
hanning_lut[i] = $rtoi(
(2.0**(DATA_WIDTH-1)-1) * 0.5 * (1.0 - $cos(2.0*3.14159265*i/(N-1)))
);
// Hamming: 0.54 - 0.46 * cos(2πn/(N-1))
hamming_lut[i] = $rtoi(
(2.0**(DATA_WIDTH-1)-1) * (0.54 - 0.46*$cos(2.0*3.14159265*i/(N-1)))
);
// Blackman: 0.42 - 0.5*cos(...) + 0.08*cos(...)
blackman_lut[i] = $rtoi(
(2.0**(DATA_WIDTH-1)-1) * (0.42 - 0.5*$cos(2.0*3.14159265*i/(N-1))
+ 0.08*$cos(4.0*3.14159265*i/(N-1)))
);
end
end
always @(*) begin
case (window_type)
2'b00: window_out = {1'b0, {(DATA_WIDTH-1){1'b1}}}; // 矩形窗=1
2'b01: window_out = hanning_lut[index];
2'b10: window_out = hamming_lut[index];
2'b11: window_out = blackman_lut[index];
endcase
end
endmodule//=============================================
// dsp_peak_detector.v
// 频谱峰值检测器
// 检测DFT输出中的频谱峰值位置和幅度
//=============================================
module dsp_peak_detector #(
parameter N = 16,
parameter DATA_WIDTH = 32
)(
input wire clk,
input wire rst_n,
input wire start,
input wire [DATA_WIDTH-1:0] spectrum [0:N-1], // 幅度谱输入
output reg [$clog2(N)-1:0] peak_index, // 峰值频率索引
output reg [DATA_WIDTH-1:0] peak_value, // 峰值幅度
output reg [$clog2(N)-1:0] second_peak_index, // 第二峰值索引
output reg [DATA_WIDTH-1:0] second_peak_value, // 第二峰值幅度
output reg done
);
localparam IDLE = 0, SCAN = 1, DONE = 2;
reg [1:0] state;
reg [$clog2(N)-1:0] scan_idx;
always @(posedge clk or negedge rst_n) begin
if (!rst_n) begin
state <= IDLE;
peak_index <= 0; peak_value <= 0;
second_peak_index <= 0; second_peak_value <= 0;
done <= 1'b0;
end else begin
case (state)
IDLE: begin
done <= 1'b0;
if (start) begin
peak_index <= 0; peak_value <= 0;
second_peak_index <= 0; second_peak_value <= 0;
scan_idx <= 0;
state <= SCAN;
end
end
SCAN: begin
if (spectrum[scan_idx] > peak_value) begin
second_peak_value <= peak_value;
second_peak_index <= peak_index;
peak_value <= spectrum[scan_idx];
peak_index <= scan_idx;
end else if (spectrum[scan_idx] > second_peak_value) begin
second_peak_value <= spectrum[scan_idx];
second_peak_index <= scan_idx;
end
if (scan_idx == N - 1)
state <= DONE;
else
scan_idx <= scan_idx + 1'b1;
end
DONE: begin
done <= 1'b1;
state <= IDLE;
end
endcase
end
end
endmodule频率分辨率 Δf = fs/N，意味着：
计算4点DFT：x[n] = {1, 2, 3, 4}，求X[k]（k=0,1,2,3）
对x[n] = cos(2π·3.5/16·n)，n=0..15，分别用矩形窗和Hanning窗做16点DFT，比较频谱泄漏情况。
修改 dsp_dft_engine.v，实现全流水线结构，使得每个时钟周期都能输出一个频点的计算结果。
使用峰值检测器的结果，实现频率插值算法（如抛物线插值），提高频率估计精度到子 bin 级别。
✅ 掌握了DFT的定义与计算方法
✅ 理解了DFT性质与循环卷积
✅ 实现了DFT计算引擎（Verilog）
✅ 实现了窗函数生成器与峰值检测
✅ 理解了频谱泄漏与分辨率的关系
离散傅里叶变换(DFT)在进阶应用中还有以下重要主题值得深入研究：
本课内容在整个DSP课程体系中处于承上启下的位置：
在实际工程中，离散傅里叶变换(DFT)的设计经验总结如下：（1）始终从系统级需求出发，不要过早陷入实现细节；（2）先用浮点仿真验证算法正确性，再转为定点；（3）系数量化后必须重新验证频率响应；（4）硬件实现前用C/MATLAB模型作为参考；（5）综合后检查时序和资源是否符合预期；（6）板级验证时用已知的测试信号对比参考输出。每一步验证都是必要的，跳过任何一步都可能导致最终的调试时间成倍增加。
在实际工程中，离散傅里叶变换(DFT)的设计经验总结如下：（1）始终从系统级需求出发，不要过早陷入实现细节；（2）先用浮点仿真验证算法正确性，再转为定点；（3）系数量化后必须重新验证频率响应；（4）硬件实现前用C/MATLAB模型作为参考；（5）综合后检查时序和资源是否符合预期；（6）板级验证时用已知的测试信号对比参考输出。每一步验证都是必要的，跳过任何一步都可能导致最终的调试时间成倍增加。
在实际工程中，离散傅里叶变换(DFT)的设计经验总结如下：（1）始终从系统级需求出发，不要过早陷入实现细节；（2）先用浮点仿真验证算法正确性，再转为定点；（3）系数量化后必须重新验证频率响应；（4）硬件实现前用C/MATLAB模型作为参考；（5）综合后检查时序和资源是否符合预期；（6）板级验证时用已知的测试信号对比参考输出。每一步验证都是必要的，跳过任何一步都可能导致最终的调试时间成倍增加。