給大家安利一個 AI 學習神站!在這個 AI 卷成紅海的時代,甭管你是硬核開發者還是代碼小白,啃透 AI?技能樹都是剛需。這站牛逼之處在于:全程用 "變量名式" 幽默 + 生活化類比拆解 AI,從入門到入土(啊不,到精通)一站式通關,理解效率直接拉滿。話不多說,鏈接甩這了→前言 – 人工智能教程
項目目標:?在FPGA上實現一個使用漢明窗設計的FIR低通濾波器,并通過仿真驗證其濾波效果。
一、FIR 原理
1. FIR = 加權移動平均(升級版): 想象在聽一串連續的聲音采樣值(x[n], x[n-1], x[n-2], ...)。FIR 濾波器計算當前輸出 y[n] 的方式是:把當前的和過去若干個輸入樣本,各自乘上一個特定的“權重”(系數 h[k]),然后把所有這些乘積加起來。
2. “有限”的含義: 它只使用有限個 (N+1 個) 過去的樣本(x[n] 到 x[n-N])。N+1 就是濾波器的抽頭數(Taps),N 是階數。抽頭數越多,濾波器能實現的頻率響應特性(比如截止陡峭度、阻帶衰減)通常越好,但計算量也越大。
3. 公式: y[n] = h[0]*x[n] + h[1]*x[n-1] + h[2]*x[n-2] + ... + h[N]*x[n-N]
4. 關鍵點: 這些系數 h[0], h[1], ..., h[N] 決定了濾波器的特性(低通、高通、帶通等)。系數不同,濾波效果就完全不同。
5. 窗函數(漢明窗)的作用: 直接計算理想濾波器的系數(無限長)不可行。窗函數(如漢明窗)像一把平滑的“剪刀”,把理想的、無限長的系數序列“截短”成我們需要的有限長 (N+1個),同時盡量減小截短帶來的不良影響(如吉布斯效應),讓濾波器的實際頻率響應更接近理想狀態。
二、MATLAB 部分
目標: 根據想要的濾波器特性(低通,Fs(采樣頻率)=100kHz, Fc(截止頻率)=10kHz),計算出 32 個 (N+1=32) 最優的量化系數 h_quant[0:31]。
步驟:
1. 計算歸一化截止頻率: Wn = Fc / (Fs/2) = 10000 / 50000 = 0.2
2. 調用 fir1 函數設置濾波器: h = fir1(31, 0.2, 'low', hamming(32));
?? ??? ?31: 濾波器階數 N (抽頭數 Taps = N+1 = 32)。
?? ??? ?0.2: 歸一化截止頻率 Wn。
?? ??? ?'low': 低通濾波器。
?? ??? ?hamming(32): 使用 32 點的漢明窗。
% 濾波器參數
Fs = 100000; % 采樣頻率 (Hz)
Fc = 10000; % 截止頻率 (Hz)
N = 31; % 濾波器階數 (抽頭數 = N + 1 = 32)
% 計算歸一化截止頻率 (范圍 0-1, 1對應 Fs/2)
Wn = Fc / (Fs/2); % Wn = 10000 / 50000 = 0.2
% 使用 fir1 函數設計漢明窗低通FIR濾波器
% fir1(階數, 歸一化截止頻率, 濾波器類型, 窗函數)
h = fir1(N, Wn, 'low', hamming(N+1)); % hamming(N+1) 生成 N+1 點的漢明窗
????????同時需要檢查圖形,通帶是否在0-10kHz相對平坦?阻帶衰減是否足夠(漢明窗典型約-53dB)?過渡帶是否在預期位置?
% 繪制頻率響應
freqz(h, 1, 1024, Fs); % freqz(系數, 分母(1表示FIR), FFT點數, 采樣頻率)
title('FIR Lowpass Filter Frequency Response (Hamming Window, N=31)');
grid on;
3. 量化系數 (關鍵!FPGA 需要整數、定點數):
?? ??? ?h 是 MATLAB 算出的浮點數系數 (范圍大約 -1 到 1)。
?? ??? ?確定系數位寬 COEFF_WIDTH (例如 16 bit)。
?? ??? ?確定量化格式: Q 格式 (常用 Q1.15,即 1 位符號位 + 15 位小數位)。
COEFF_WIDTH = 16;
FRAC_BITS = COEFF_WIDTH - 1; % 15 for Q1.15
h_quant = round(h * (2^FRAC_BITS)); % 放大 2^15 倍后四舍五入取整
% 飽和處理 (防止溢出)
max_val = 2^(COEFF_WIDTH-1) - 1; % 32767 for 16-bit
min_val = -2^(COEFF_WIDTH-1); % -32768 for 16-bit
h_quant(h_quant > max_val) = max_val;
h_quant(h_quant < min_val) = min_val;
4. 生成 .coe 文件 (供 Vivado ROM IP 使用):
% 打開文件用于寫入
fid = fopen('fir_coefficients.coe', 'w');
% 寫入.coe文件頭 (指定Radix和系數值)
fprintf(fid, 'memory_initialization_radix = 10;\n'); % 或 2, 16
fprintf(fid, 'memory_initialization_vector = \n');
% 寫入系數 (十進制整數形式,最后一個系數后面不加逗號)
for i = 1:length(h_quant)if i < length(h_quant)fprintf(fid, '%d,\n', h_quant(i));elsefprintf(fid, '%d;\n', h_quant(i)); % 最后一個系數以分號結尾end
end
fclose(fid);
disp('COE file "fir_coefficients.coe" generated.');
三、vivado?工程
整體框圖如下:
1. 創建用于存儲系數的ROM (Block Memory Generator IP)
2. 實現串行FIR濾波器結構
核心思想為:
-
使用單個乘法器依次計算每個抽頭的乘積
-
通過狀態機精確控制計算流程
-
移位寄存器存儲歷史樣本
-
ROM存儲濾波器系數
完整工作流程為:
-
復位階段:所有寄存器清零,狀態機歸零
-
新樣本到達:進入狀態0
-
移位寄存器更新
-
累加器和地址初始化
-
-
乘累加階段:狀態1到TAPS-1
-
依次計算每個抽頭的乘積
-
累加前一個抽頭的計算結果
-
-
輸出階段:狀態TAPS
-
輸出累加結果
-
返回狀態0等待新樣本
-
`timescale 1ns / 1psmodule serial_fir#(parameter DATA_WIDTH = 12,parameter OUTPUT_WIDTH = 28,parameter COEFF_WIDTH = 16,parameter TAPS = 32
)(input wire clk,input wire rst,input wire signed [DATA_WIDTH - 1 : 0] data_in,output reg signed [OUTPUT_WIDTH - 1 : 0] data_out
);
//移位寄存器鏈;存儲當前與過去的輸入樣本reg signed [DATA_WIDTH-1:0] shift_reg [0:TAPS-1];//前面是每個元素的位寬,后面是一個數組及索引范圍integer i;//FIR系數接收wire signed[COEFF_WIDTH-1:0] coeff_out;reg [$clog2(TAPS)-1:0] addr ;fir_coeff_rom your_instance_name (.clka(clk), // input wire clka.addra(addr), // input wire [4 : 0] addra.douta(coeff_out) // output wire [15 : 0] douta
);//乘法器(為了高時鐘F加流水線級),累加器,狀態寄存器reg signed[OUTPUT_WIDTH-1:0] accumulator;reg signed[$clog2(TAPS):0] state;reg signed[OUTPUT_WIDTH-1:0]product_reg;reg signed[OUTPUT_WIDTH-1:0]mult_result;always@(posedge clk or posedge rst)if(rst)beginfor(i=0; i<TAPS; i=i+1) shift_reg[i] <= 0;addr <= 0;accumulator <= 0;state <= 0;product_reg <= 0;data_out <=0;end else begincase(state)0: beginfor(i=TAPS-1; i>0; i=i-1)beginshift_reg[i] <= shift_reg[i-1];endshift_reg[0] <= data_in;state <= 1;enddefault beginif(state <= TAPS) beginmult_result <= shift_reg[addr] * coeff_out;product_reg <= mult_result;accumulator <= accumulator + product_reg;addr <= addr + 1;state <= state + 1;end else begindata_out <= accumulator;state <= 0;endendendcaseendendmodule
四、使用乘法器IP核進行硬件加速
????????雖然Verilog的 * 操作符在小位寬時綜合工具可能能推斷出乘法器,但對于16位或更大位寬,使用Xilinx專用的DSP Slice IP核 (Multiplier) 通常性能更好、資源利用更優、時序更容易滿足。
修改 serial_fir.v 代碼:
? ? ? ? 1. 找到原來的乘法操作 shift_reg[addr] * coeff_out;
? ? ? ? 2. 同時將原來的代碼進行時序狀態機優化,
? ? ? ? 2. 替換為乘法器IP核實例化:
`timescale 1ns / 1psmodule serial_fir#(parameter DATA_WIDTH = 12,parameter OUTPUT_WIDTH = 28,parameter COEFF_WIDTH = 16,parameter TAPS = 32
)(input wire clk,input wire rst,input wire signed [DATA_WIDTH - 1 : 0] data_in,output reg signed [OUTPUT_WIDTH - 1 : 0] data_out
);// 移位寄存器鏈reg signed [DATA_WIDTH-1:0] shift_reg [0:TAPS-1];integer i;// FIR系數接收wire signed [COEFF_WIDTH-1:0] coeff_out;reg [$clog2(TAPS)-1:0] addr;// 狀態寄存器 - 需要覆蓋0到TAPS+1狀態reg [$clog2(TAPS+2):0] state; // 修正位寬// 乘法器、累加器相關reg signed [OUTPUT_WIDTH-1:0] accumulator;reg signed [OUTPUT_WIDTH-1:0] product_reg;wire signed [OUTPUT_WIDTH-1:0] mult_result; // 改為wire類型// 系數ROMfir_coeff_rom coeff_rom_inst (.clka(clk),.addra(addr),.douta(coeff_out));// 乘法器IP核實例化signed_multiplier mult_inst (.CLK(clk),.A(shift_reg[addr]),.B(coeff_out),.P(mult_result));// 狀態定義localparam S_LOAD = 0; // 加載新樣本localparam S_MULT = 1; // 啟動乘法localparam S_ACC = 2; // 累加結果localparam S_OUTPUT = 3; // 輸出結果always @(posedge clk or posedge rst) beginif (rst) begin// 復位所有寄存器for (i = 0; i < TAPS; i = i + 1) shift_reg[i] <= 0;addr <= 0;accumulator <= 0;state <= S_LOAD;product_reg <= 0;data_out <= 0;end else begincase (state)S_LOAD: begin// 移位寄存器更新for (i = TAPS-1; i > 0; i = i - 1) beginshift_reg[i] <= shift_reg[i-1];endshift_reg[0] <= data_in;// 初始化累加器和地址accumulator <= 0;addr <= 0;// 進入乘法狀態state <= S_MULT;endS_MULT: begin// 啟動乘法(結果將在下一周期出現在mult_result)// 不需要直接賦值,乘法器IP會自動計算// 存儲上一個乘法結果(如果有)// 對于第一個抽頭,product_reg為0// 進入累加狀態state <= S_ACC;endS_ACC: begin// 累加前一個乘法結果accumulator <= accumulator + product_reg;// 存儲當前乘法結果(用于下一個抽頭的累加)product_reg <= mult_result;// 更新地址addr <= addr + 1;// 判斷是否完成所有抽頭if (addr < TAPS - 1) beginstate <= S_MULT; // 繼續下一個抽頭end else beginstate <= S_OUTPUT; // 所有抽頭處理完成endendS_OUTPUT: begin// 累加最后一個抽頭的結果accumulator <= accumulator + product_reg;// 輸出最終結果data_out <= accumulator;// 返回加載狀態state <= S_LOAD;endendcaseendend
endmodule
五、編寫仿真文件查看波形
`timescale 1ns/1psmodule serial_fir_tb();localparam DATA_WIDTH = 12;
localparam OUTPUT_WIDTH = 28;
localparam COEFF_WIDTH = 16;
localparam TAPS = 32;
localparam CLK_PERIOD = 10;reg clk;
reg rst;
reg signed [DATA_WIDTH - 1 : 0] data_in;
wire signed [OUTPUT_WIDTH - 1 : 0] data_out ;parameter FS_REAL = 100000;
parameter REAL_PERIOD = 1.0e9/FS_REAL;
integer t;
real time_val;
parameter freq1 = 5000;
parameter freq2 = 30000;
parameter pi = 3.1415926;
parameter amplitude = (1 << (DATA_WIDTH-1))-1;//時鐘生成
always #(CLK_PERIOD/2) clk = ~clk;//實例化模塊
serial_fir#(.DATA_WIDTH(DATA_WIDTH),.OUTPUT_WIDTH(OUTPUT_WIDTH),.COEFF_WIDTH(COEFF_WIDTH),.TAPS(TAPS)
) uut(.clk(clk),.rst(rst),.data_in(data_in),.data_out(data_out)
);//初始化
initial beginclk = 0;rst = 1;data_in = 0;#(CLK_PERIOD*2) rst = 0;//測試1:單位沖擊響應(驗證濾波系數)
data_in = (1 << (DATA_WIDTH-2));
#(CLK_PERIOD);
data_in = 0;
#(CLK_PERIOD*(TAPS*10));//測試2:混合頻率信號(驗證濾波效果)
//localparam FS_REAL = 100000.0;
//localparam REAL_PERIOD = 1.0e9 / FS_REAL;
//integer t;
//real time_val;
//real freq1 = 5000;
//real freq2 = 30000;
//real pi = 3.1415926;
//real amplitude = (1 << (DATA_WIDTH-1))-1;for(t=0; t<500; t=t+1) begintime_val = t*(REAL_PERIOD/1.0e9);//第t個樣本對應的實際時間sdata_in = $rtoi(amplitude * (0.6 * $sin(2 * pi * freq1 * time_val) + 0.4 * $sin(2 * pi * freq2 * time_val)));$display("t = %d, time_val = %f, data_in = %d, data_out = %d", t, time_val, data_in, data_out);#(REAL_PERIOD);
end//結束仿真
#(CLK_PERIOD*100);
$finish;end//波形記錄
initial begin$dumpfile("waveform.vcd");$dumpvars(0,uut);
endendmodule
仿真結果如下圖所示:
六、MATLAB進行FFT分析
將生成的dumpfile文件復制到txt文檔然后導入MATLAB進行FFT分析(一部分如下)
MATLAB代碼如下:
% FIR濾波器頻譜分析完整代碼
% 功能:讀取數據文件,計算并繪制data_in和data_out的頻譜,驗證濾波效果% --------------------------
% 1. 讀取數據文件
% --------------------------
% 確保數據文件與當前MATLAB腳本在同一目錄,或使用絕對路徑
data_filename = 'fir_data_clean.txt'; % 數據文件名
data = load(data_filename); % 讀取空格分隔的純數值文件% 提取各列數據
t = data(:, 1); % 時間索引
data_in = data(:, 2); % 輸入信號
data_out = data(:, 3); % 輸出信號% 驗證數據讀取結果
fprintf('成功讀取數據,共 %d 個樣本\n', length(t));
fprintf('時間范圍: t = %d 至 %d\n', t(1), t(end));% --------------------------
% 2. 配置信號參數
% --------------------------
FS = 100000; % 采樣率 (Hz),與仿真中的FS_REAL一致
N = length(data_in); % 數據點數
Nyquist = FS / 2; % 奈奎斯特頻率 (Hz)% --------------------------
% 3. 計算FFT并處理頻譜
% --------------------------
% 對輸入信號計算FFT
dc_in = mean(data_in); % 計算直流分量
fft_in = fft(data_in - dc_in); % 去除直流分量后計算FFT
fft_mag_in = abs(fft_in) / N * 2; % 幅度歸一化(雙邊譜轉單邊譜)
fft_mag_in(1) = fft_mag_in(1) / 2; % 修正直流分量幅度% 對輸出信號計算FFT
dc_out = mean(data_out);
fft_out = fft(data_out - dc_out);
fft_mag_out = abs(fft_out) / N * 2;
fft_mag_out(1) = fft_mag_out(1) / 2;% 生成頻率軸(單位:Hz)
freq = (0:N-1) * FS / N;% --------------------------
% 4. 繪制時域波形
% --------------------------
figure('Name', '時域波形', 'Position', [100, 100, 1000, 600]);
subplot(2, 1, 1);
plot(t, data_in, 'b-', 'LineWidth', 1);
xlabel('時間索引 t');
ylabel('幅度');
title('輸入信號 data_in 時域波形');
grid on;subplot(2, 1, 2);
plot(t, data_out, 'r-', 'LineWidth', 1);
xlabel('時間索引 t');
ylabel('幅度');
title('輸出信號 data_out 時域波形');
grid on;
sgtitle('時域波形對比');% --------------------------
% 5. 繪制頻譜圖(0至奈奎斯特頻率)
% --------------------------
figure('Name', '頻譜對比', 'Position', [200, 200, 1000, 600]);% 輸入信號頻譜
subplot(2, 1, 1);
plot(freq(1:N/2), fft_mag_in(1:N/2), 'b-', 'LineWidth', 1.2);
xlabel('頻率 (Hz)');
ylabel('幅度');
title('data_in 頻譜(5kHz + 30kHz 混合信號)');
grid on;
xlim([0, Nyquist]); % 顯示0至50kHz
hold on;
% 標記關鍵頻率
plot([5000, 5000], ylim, 'r--', 'LineWidth', 1); % 5kHz(目標保留頻率)
plot([30000, 30000], ylim, 'g--', 'LineWidth', 1); % 30kHz(目標衰減頻率)
legend('信號幅度', '5kHz(通帶)', '30kHz(阻帶)', 'Location', 'best');% 輸出信號頻譜
subplot(2, 1, 2);
plot(freq(1:N/2), fft_mag_out(1:N/2), 'r-', 'LineWidth', 1.2);
xlabel('頻率 (Hz)');
ylabel('幅度');
title('data_out 頻譜(濾波后)');
grid on;
xlim([0, Nyquist]);
hold on;
plot([5000, 5000], ylim, 'r--', 'LineWidth', 1);
plot([30000, 30000], ylim, 'g--', 'LineWidth', 1);
legend('信號幅度', '5kHz(通帶)', '30kHz(阻帶)', 'Location', 'best');sgtitle('輸入輸出信號頻譜對比');% --------------------------
% 6. 定量分析關鍵頻率成分
% --------------------------
% 查找5kHz附近的頻率索引
[idx_5k, ~] = min(abs(freq - 5000));
mag_in_5k = fft_mag_in(idx_5k); % 輸入信號5kHz幅度
mag_out_5k = fft_mag_out(idx_5k); % 輸出信號5kHz幅度% 查找30kHz附近的頻率索引
[idx_30k, ~] = min(abs(freq - 30000));
mag_in_30k = fft_mag_in(idx_30k); % 輸入信號30kHz幅度
mag_out_30k = fft_mag_out(idx_30k);% 輸出信號30kHz幅度% 計算衰減比和增益
attenuation_ratio = mag_in_30k / mag_out_30k; % 30kHz衰減倍數(越大越好)
gain_5k = mag_out_5k / mag_in_5k; % 5kHz增益(應接近1)% 打印分析結果
fprintf('\n===== 頻譜定量分析結果 =====\n');
fprintf('5kHz信號幅度: 輸入=%.2f, 輸出=%.2f\n', mag_in_5k, mag_out_5k);
fprintf('30kHz信號幅度: 輸入=%.2f, 輸出=%.2f\n', mag_in_30k, mag_out_30k);
fprintf('30kHz衰減比: %.2f 倍(理想低通濾波器應遠大于1)\n', attenuation_ratio);
fprintf('5kHz增益: %.2f(理想應接近1)\n', gain_5k);% 判斷濾波效果
if attenuation_ratio > 5 && gain_5k > 0.5fprintf('\n結論: 濾波效果符合預期,30kHz高頻被有效衰減,5kHz低頻保留良好。\n');
elsefprintf('\n結論: 濾波效果不理想,需檢查FIR濾波器系數或內部邏輯。\n');
end
從結果來看,濾波效果符合預期,原因如下:
1. 頻譜分析
- 輸入頻譜:
data_in
?在 5kHz 和 30kHz 處均有明顯峰值,符合 “5kHz + 30kHz 混合信號” 的設計預期。 - 輸出頻譜:
data_out
?僅保留了 5kHz 處的峰值,30kHz 處的峰值幾乎完全被衰減,說明 FIR 濾波器成功實現了 “保留 5kHz 低頻、衰減 30kHz 高頻” 的低通濾波功能。
2. 結論
當前頻譜結果表明,FIR 濾波器的濾波效果符合設計目標,30kHz 高頻被有效抑制,5kHz 低頻被完整保留。
后續改進方向:可向并行設計靠攏