本文借鑒 小墨學FPGA 的文章 《零基礎學FPGA（ 三十二） 寫在京城，多級FIR半帶濾波器的FPGA實現 》
地址是： http://bbs.elecfans.com/forum.php?mod=viewthread&tid=509022
假設我們的合成信號是由100,300和600hz的正弦波相加而成的信號，並將其進行8bit量化，初始采樣頻率為6400hz，現在我需要將采樣信號降為800hz，並要求獲得100hz有用信號，要求濾波器的通阻帶容限為0.1，設計這樣一個FIR半帶濾波器。
clear all
clc
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 通帶截止頻率近似3KHz，阻帶截止頻率近似5KHz
fp = 100; % 要求最后一級濾波器的通帶截止頻率
fc = 300; % 要求最后一級濾波器的阻帶截止頻率
fs = 6400; %原始采樣頻率
fs_out = 800; % 要求降頻后的采樣頻率
dev=0.1; %要求的濾波器通阻帶容限
D=fs/fs_out; %求總抽取倍數
N=log2(D) ; %求濾波器級數
devi=dev/N; % 求每級濾波器通阻帶容限
% 設計前2級半帶濾波器
hin1 = firhalfband('minorder',fp/(fs/2),dev); %designs a lowpass minimum-order filter,
% with passband edge fp. The peak ripple is constrained by the scalar dev.
% This design uses the equiripple method.
hin2 = firhalfband('minorder',fp/(fs/4),dev);
fid=fopen('E:\WORK\DRS6000_Q\MATLAB CODE\hin1.txt','wt');
fprintf(fid, '%d\n', hin1);
fclose(fid);
fid=fopen('E:\WORK\DRS6000_Q\MATLAB CODE\hin2.txt','wt');
fprintf(fid, '%d\n', hin2);
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 設計最后一級最優低通濾波器
f = [fp,fc]; % 最優低通濾波器帶寬
a = [1 0]; % 最優低通濾波器幅頻特性，低通濾波器
devk = [devi,dev] ; % 要求的濾波器通阻帶容限
mag = [1 1 0 0]; %最優濾波器幅度特性
fpm = [0 fp/((fs_out*2)/2) fc/((fs_out*2)/2) 1];
[n,wn,beta,firtype] = kaiserord(f,a,devk,fs_out*2); %調用函數求濾波器最小階數
hin3 =firpm(n,fpm,mag); %調用函數設計最優濾波器
% freqz(hin3);
fid=fopen('E:\WORK\DRS6000_Q\MATLAB CODE\hin3.txt','wt');
fprintf(fid, '%d\n', hin3);
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f1 = 100; % 信號頻率
f2 = 300; % 信號頻率
f3 = 600; % 信號頻率
lengthy = 1024;
n = 0:1:lengthy-1;
t = n/fs;
si = sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t);
figure
plot(t,si);
freq=linspace(-fs/2,fs/2,lengthy);
figure
plot(freq,fftshift( abs(fft(si)) ));
xlim( [ 0, 1000 ] ); %設置橫軸范圍0~1000
xlabel('Frequency ');
title('Magnitude spectrum of chirp signal');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%濾波與抽取過程
F1_out = filter (hin1,1,si); %第一級半帶濾波器
s1 = F1_out (1:2:length(F1_out)); %2倍抽取
F1_max = max (abs(F1_out));
F2_out = filter (hin2,1,s1); %第二級半帶濾波器
s2 = F2_out (1:2:length(F2_out)); %2倍抽取
F3_out = filter (hin3,1,s2); %第二級半帶濾波器
sout = F3_out (1:2:length(F3_out)); %2倍抽取
sout = sout/max(sout);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lengths = lengthy/8;
n1 = 0:1:lengths-1;
t1=n1/fs;
figure
plot(t1,sout);
freq=linspace(-fs/16,fs/16,lengths);
figure
plot(freq,fftshift( abs(fft(sout)) ));
xlim( [ 0, 1000 ] ); %設置橫軸范圍0~1000
xlabel('Frequency ');
title('Magnitude spectrum of chirp signal');
figure
subplot(211)
stem(si);
subplot(212)
stem(sout);
% % 前2級半帶濾波器的幅頻相應
% figure
% impz(hin1);
% figure
% impz(hin2);
%
% h1 = mfilt.firdecim(2,hin1); % Create a polyphase decimator
% h2 = mfilt.firdecim(2,hin2); % Create a polyphase decimator
% freqz(h1);
% freqz(h2);