NB. Fourier regression NB. DeLurgio P250 NB. * x reg_fourier y NB. * x estim_fourier y NB. plot |: x estim_fourier y NB. ============================= require 'plot numeric trig' NB. ============================== reg_fourier=: 4 : 0 NB. x is jisuu JISUU=. x NB. usage: jisuu(e.g. 5) reg_fourier y //e.g. DATF 'NR Y0'=. take_12time y TREND0=.(find_trend0 Y0) Y1=. ({. TREND0)+ ({: TREND0)* >: i. # Y0 (Y0-Y1) %. (+: JISUU){."1 (NR;JISUU) reg_fourier_sub Y0 NB. sample -->12{. "1 occure error ) take_12time=: 3 : 0 NR=. <. (#y) % 12 Y0=. (NR * 12){. y NR;Y0 ) reg_fourier_sub=: 4 : 0 'NR JISUU'=: x NB. x is NR T=. >: i.# y W=. 2p1 * NR %# Y0 NB. 48 or 12 times TMPA=. ,.(T*W)* L:0 {@> >: i. >: JISUU TMPS0=.(cos L:0 TMPA), sin L:0 TMPA NR=. # TMPS0 IND=.,|:((NR %-: NR),-: NR)$ i. NR TMPA1=. |: ;("1) IND{TMPS0 ) find_trend0=: 3 : 0 B0=.12%~(+/%#) -/"1 (12}. y),. _12}. y NB. bo A0=.((+/%#) y)-B0*(-:>:#y) REG=:A0,B0 ) estim_fourier=: 4 : 0 NB. plot |: x estim_fourier DATF JISUU=: x 'NR Y0'=. take_12time y T=. >: i.# Y0 W=. 2p1 * NR % # Y0 NB. 48 or 12 times TREND0=. find_trend0 Y0 FOURIER=: JISUU reg_fourier Y0 NB. ---------------------- TMPREG=. ({. TREND0)+({: TREND0)* T FOURIER1=.(,2 ,JISUU) $ (,|: (,JISUU,2) $ i.+: JISUU){ FOURIER NB. limit 10 TMPCOS=.+/ ({. FOURIER1)*/ cos 2p1 *(NR*T)%(#T) TMPSIN=.+/ ({: FOURIER1)*/ sin 2p1 *(NR*T)%(#T) Y0,.TMPREG,. TMPREG + TMPCOS + TMPSIN NB. actual trend fourier ) NB. =================Caliculas========= NB. Fourier series fou=: 2&o.@(*& 0 1 2)"0 mp=: +/ . * m=: i. 3 3 L=: m&mp LO=: L@ NB. from lab of FFT NB. Discrete Fourier translation dfft=: 3 : '+/ y * ^ (#y) %~ (- o. 0j2 ) * */~ i.#y' NB. 'bar' plot dfft y is DATF