%《數字信號處理教程——MATLAB釋義與實現》第七章例7.6.2程序hc762
% 希爾波特變換器的設計
% 電子工業出版社出版  陳懷琛編著 2004年9月
%
N=input('濾波器長度＝(輸入奇數） '); % 給出濾波器長度，它等于階數加一
tau=(N-1)/2;n=[0:2*tau]+1e-10;              % 給出n數組，加微小偏移
hd=2*(sin((n-tau)*pi/2)).^2./(n-tau)/pi;    % 矩形窗截斷的系數向量
hhn=hd.*hann(N)';                  % 漢寧窗截斷的系數向量
[H,w]=freqz(hd,1,'whole');              % 求濾波器全頻域頻率響應
[Hhn,w]=freqz(hhn,1,'whole');
subplot(2,2,1),stem([0:N-1],hd,'.')     % 畫出系數向量hd
axis([0,N-1,-1,1]),grid on
subplot(2,2,2),plot(w,abs(H)),grid on   % 畫出幅頻特性|H|
subplot(2,2,3),plot(w,unwrap(angle(H)))  % 畫出相頻頻特性
axis([2,4,-2*(N-1),-(N-1)])
hold on,plot(w,-tau*w,':'),grid on
subplot(2,2,4),plot(w,abs(Hhn)),grid on % 畫出幅頻特性|Hhn|
set(gcf,'color','w'),pause              % 置圖形背景色為白
% 用信號處理工具箱中的最優化設計函數設計希爾波特變換器
F=[0,1],A=[1,1]; 
b=firls(N-1,F,A,'Hilbert')          % 進行最小二乘設計
b1=remez(N-1,F,A,'Hilbert')         % 進行最優等波動設計