
function h = lpdiff(K,N)

% h = lpdiff(K,N)
% lowpass differentiator
% N : length of impulse response h
% K : flatness (approximate) at DC (1 <= K <= 8)

a = zeros(1,K);

if K == 1
    a(1) = 1/2;
elseif K == 2
    a(1) = 1/9*pi^2-1/6;
    a(2) = 2/9*pi^2-4/3;
elseif K == 3
    a(1) = -13/288*pi^2+1/48+7/480*pi^4;
    a(2) = -16/9*pi^2+16/15+14/75*pi^4;
    a(3) = -81/32*pi^2+243/80+189/800*pi^4;
elseif K == 4
    a(1) = -1/720+29/4320*pi^2-427/64800*pi^4+31/18900*pi^6;
    a(2) = -16/45+208/135*pi^2-2366/2025*pi^4+496/4725*pi^6;
    a(3) = -2187/560+2187/160*pi^2-5103/800*pi^4+2511/4900*pi^6;
    a(4) = -2048/315+2048/135*pi^2-12544/2025*pi^4+15872/33075*pi^6;
elseif K == 5
    a(2) = 1.996550426365710869302158146;
    a(4) = 2.7661556746480243935701781;
    a(3) = 2.85469852315917092007298400;
    a(5) = 1.1912526731983849661041972;
    a(1) = .999997422818794233131789510849;
elseif K == 6
    a(1) = .999999965390251190899080316192;
    a(3) = 2.9798121081274202540352886;
    a(2) = 1.99980491783061681143438898;
    a(6) = 1.164693736431702890703629;
    a(5) = 3.072298081103930132715621;
    a(4) = 3.654848810665006991811345;
elseif K == 7
    a(6) = 3.25155058776190389487903206011674797453005622338256583;
    a(7) = 1.10579058782215448021144271382362003105778593441657857;
    a(2) = 1.99999209478894973314698867639078680064709305491874928549;
    a(1) = .999999999663375885647929058808192901256576143562055933204249;
    a(4) = 3.933640382058333295378017583093438408874386149124063580;
    a(3) = 2.998026251788322276959914483777903388049932800791857886;
    a(5) = 4.34009536650982682211770139544385624484028921852173053;
elseif K == 8
    a(2) = 1.9999997592548908941351745;
    a(8) = 1.0276795254305789375815;
    a(5) = 4.8392625101024630589032;
    a(6) = 4.8995532200454177536038407022;
    a(7) = 3.324490070500423334579496514;
    a(3) = 2.99985691877486768420045;
    a(1) = .999999999997520082999046359037;
    a(4) = 3.99071975023105221830281;
end

n = 0:N-1;
k = 1;
h = a(1)*sin(2*pi/N*k*(n-(N-1)/2));
for k = 2:K
    h = h + a(k)*sin(2*pi/N*k*(n-(N-1)/2));
end
h = h*4/N/N;

