最小二乗法によるFIRフィルタの設計

clear; close all;

N = 61;
A = 0;
b = 0;
syms w;
M = (N-1)/2;
c = [w];%dummy
for n=1:M
    c(end+1) = 2*cos(w*n);
end
c(1) = 1;
c = fliplr(c);
c = transpose(c);

omega = linspace(0, pi, N+1);
%遮断周波数の定義
delta_w = 0.1*pi;%転移域の幅
wc = 0.5*pi;%遮断周波数
wd = wc+delta_w;
%重み関数Wの作成
W = zeros(size(omega));
Wconst = 0.9;%任意の値 
W(find(omega<wc)) = 1;
W(find(omega>=wc & omega<wd)) = 0;
W(find(omega>wd))=Wconst;
%ディジタルLPFの理想特性
D = zeros(size(omega));
D(find(omega<wc))=1;
D(find(omega>wc))=0;
A = zeros(length(c),length(c));
b = zeros(size(c));

K = N;%Kは任意の値（十分に大きく取る）
for k=0:K-1
    k    
    ck_quad = subs(c,w,omega(k+1))*transpose(subs(c,w,omega(k+1)));
    A = A + eval(W(k+1)*ck_quad);
    b = b + eval(W(k+1)*D(k+1)*subs(c,w,omega(k+1)));
end
h = A\b;
Hw = transpose(c)*h;
Hw_val = eval(subs(Hw,w,omega));
Hw_abs_normal = abs(Hw_val)/max(abs(Hw_val));

figure();
plot(omega, Hw_abs_normal);
title('最小二乗法による直線位相FIRフィルタの振幅特性');

figure();
plot(omega, 20*log10(Hw_abs_normal));
title('最小二乗法による直線位相FIRフィルタ振幅特性（dB表示）');

big;

N=61の場合