* * UBICsignal2rfcg.f - gateway function for UBICsignal2rfc.f * subroutine mexfunction(nlhs, plhs, nrhs, prhs) *----------------------------------------------------------------------- implicit none integer*8 plhs(*), prhs(*) Integer*8 pscr integer*4 nlhs, nrhs * integer*8 mxcreatefull, mxgetm, mxgetn, mxgetpr integer*8 pr_xs, pr_x integer*8 pr_wcycles, pr_wmean, pr_wtime integer*8 pr_hcycles, pr_hmean, pr_htime integer*8 pr_T integer m_x, n_x, len_xs, len_x, len_w, len_h *----------------------------------------------------------------------- * check for proper number of arguments * if (nrhs .ne. 1) then call mexerrmsgtxt('UBIICsignal2rfc requires one input argument') elseif (nlhs .lt. 1) then call mexerrmsgtxt & ('UBICsignal2rfc requires at least one output argument') elseif (nlhs .gt. 7) then call mexerrmsgtxt & ('UBICsignal2rfc requires max four output arguments') endif * * dimensions of x * m_x = mxgetm(prhs(1)) n_x = mxgetn(prhs(1)) len_x=m_x*n_x * * create a matrix for calculation * pr_x = mxgetpr(prhs(1)) pr_xs = mxgetpr(mxcreatefull(len_x,1,0)) call signal_simplify(%val(pr_x),len_x,%val(pr_xs),len_xs) * * create matricies for return arguments * plhs(1) = mxcreatefull(len_xs,1,0) plhs(2) = mxcreatefull(len_xs,1,0) plhs(3) = mxcreatefull(len_xs,1,0) plhs(4) = mxcreatefull(len_xs,1,0) plhs(5) = mxcreatefull(len_xs,1,0) plhs(6) = mxcreatefull(len_xs,1,0) pscr = mxcreatefull(len_xs,1,0) * * assign pointers to the various parameters * pr_wcycles = mxgetpr(plhs(1)) pr_hcycles = mxgetpr(plhs(2)) pr_wmean = mxgetpr(plhs(3)) pr_hmean = mxgetpr(plhs(4)) pr_wtime = mxgetpr(plhs(5)) pr_htime = mxgetpr(plhs(6)) pr_T = mxgetpr(pscr) * * do the rfc * call UBICsignal2rfc(%val(pr_xs),len_xs, & %val(pr_wcycles),%val(pr_wmean),%val(pr_wtime),len_w, & %val(pr_hcycles),%val(pr_hmean),%val(pr_htime),len_h, & %val(pr_T)) * if (n_x.eq.1) then call mxsetn(plhs(1),1) call mxsetm(plhs(1),len_w) call mxsetn(plhs(2),1) call mxsetm(plhs(2),len_h) call mxsetn(plhs(3),1) call mxsetm(plhs(3),len_w) call mxsetn(plhs(4),1) call mxsetm(plhs(4),len_h) call mxsetn(plhs(5),1) call mxsetm(plhs(5),len_w) call mxsetn(plhs(6),1) call mxsetm(plhs(6),len_h) else call mxsetm(plhs(1),1) call mxsetn(plhs(1),len_w) call mxsetm(plhs(2),1) call mxsetn(plhs(2),len_h) call mxsetm(plhs(3),1) call mxsetn(plhs(3),len_w) call mxsetm(plhs(4),1) call mxsetn(plhs(4),len_h) call mxsetm(plhs(5),1) call mxsetn(plhs(5),len_w) call mxsetm(plhs(6),1) call mxsetn(plhs(6),len_h) endif C call mxfreematrix(pr_xs) C call mxfreematrix(pr_wcycles+len_w+1) C call mxfreematrix(pr_wmean+len_w+1) C call mxfreematrix(pr_hcycles+len_h+1) C call mxfreematrix(pr_hmean+len_h+1) * return end c----------------------------------------------------------------------- subroutine UBICsignal2rfc(x, len_x, wcycles, wmean, wtime, len_w, & hcycles, hmean, htime, len_h, T) implicit none integer len_x, len_w, len_h real*8 T(1) real*8 x(1) real*8 wcycles(1), wmean(1), wtime(1) real*8 hcycles(1), hmean(1), htime(1) * integer n, i, j, k, next real*8 e(1024), x0, y0 * n=1 j=0 k=0 next=1 E(1)=x(1) T(1)=1 do while (next.le.len_x) if (N.lt.3) then N=N+1 next=next+1 if (next.gt.len_x) goto 99 E(N)=x(next) T(N)=next else X0=abs(E(N)-E(N-1)) Y0=abs(E(N-1)-E(N-2)) if (X0.lt.Y0) then N=N+1 next=next+1 if (next.gt.len_x) goto 99 E(N)=x(next) T(N)=next else if (N.eq.3) then k=k+1 hcycles(k)=Y0 hmean(k)=(E(N-1)+E(N-2))*0.5 htime(k)=ABS(T(N-1)-T(N-2)) E(1)=E(2) E(2)=E(3) T(1)=T(2) T(2)=T(3) N=2 else j=j+1 wcycles(j)=Y0 wmean(j)=(E(N-1)+E(N-2))*0.5 wtime(j)=ABS(T(N-1)-T(N-2)) N=N-2 E(N)=E(N+2) T(N)=T(N+2) endif endif endif enddo 99 continue if (N.gt.1) then N=N-1 do i=1,N-1 hcycles(i+k)=abs(E(i)-E(i+1)) hmean(i+k)=(E(i)+E(i+1))*0.5 htime(i+k)=ABS(T(i)-T(i+1)) enddo endif len_w=j len_h=i+k-1 call sort2(len_w,wcycles,wmean,wtime) call sort2(len_h,hcycles,hmean,htime) end C----------------------------------------------------------------------- subroutine signal_simplify(x, len_x, xs, len_xs) implicit none integer len_x, len_xs real*8 x(1), xs(1) * integer i, up * xs(1)=x(1) i=2 do while (x(i).eq.x(i-1)) i=i+1 enddo if (x(i).gt.x(i-1)) then up=1 else up=0 endif len_xs=2 do i=2,len_x-1 if (up.eq.1) then if (x(i-1).lt.x(i).and.x(i).gt.x(i+1)) then xs(len_xs)=x(i) len_xs=len_xs+1 up=0 endif else if (x(i-1).gt.x(i).and.x(i).lt.x(i+1)) then xs(len_xs)=x(i) len_xs=len_xs+1 up=1 endif endif enddo xs(len_xs)=x(len_x) end c----------------------------------------------------------------------- subroutine sort2(n,arr,brr,crr) * sort from larger to smaller implicit none integer n real*8 arr(n),brr(n),crr(n) * integer m,nstack parameter (m=16,nstack=1024) integer i,ir,j,jstack,k,l,istack(nstack) real*8 a,b,c,temp jstack=0 l=1 ir=n 1 if(ir-l.lt.m)then do 12 j=l+1,ir a=arr(j) b=brr(j) c=crr(j) do 11 i=j-1,l,-1 if(arr(i).ge.a)goto 2 arr(i+1)=arr(i) brr(i+1)=brr(i) crr(i+1)=crr(i) 11 continue i=l-1 2 arr(i+1)=a brr(i+1)=b crr(i+1)=c 12 continue if(jstack.eq.0)return ir=istack(jstack) l=istack(jstack-1) jstack=jstack-2 else k=(l+ir)/2 temp=arr(k) arr(k)=arr(l+1) arr(l+1)=temp temp=brr(k) brr(k)=brr(l+1) brr(l+1)=temp temp=crr(k) crr(k)=crr(l+1) crr(l+1)=temp if(arr(l).lt.arr(ir))then temp=arr(l) arr(l)=arr(ir) arr(ir)=temp temp=brr(l) brr(l)=brr(ir) brr(ir)=temp temp=crr(l) crr(l)=crr(ir) crr(ir)=temp endif if(arr(l+1).lt.arr(ir))then temp=arr(l+1) arr(l+1)=arr(ir) arr(ir)=temp temp=brr(l+1) brr(l+1)=brr(ir) brr(ir)=temp temp=crr(l+1) crr(l+1)=crr(ir) crr(ir)=temp endif if(arr(l).lt.arr(l+1))then temp=arr(l) arr(l)=arr(l+1) arr(l+1)=temp temp=brr(l) brr(l)=brr(l+1) brr(l+1)=temp temp=crr(l) crr(l)=crr(l+1) crr(l+1)=temp endif i=l+1 j=ir a=arr(l+1) b=brr(l+1) c=crr(l+1) 3 continue i=i+1 if(arr(i).gt.a)goto 3 4 continue j=j-1 if(arr(j).lt.a)goto 4 if(j.lt.i)goto 5 temp=arr(i) arr(i)=arr(j) arr(j)=temp temp=brr(i) brr(i)=brr(j) brr(j)=temp temp=crr(i) crr(i)=crr(j) crr(j)=temp goto 3 5 arr(l+1)=arr(j) arr(j)=a brr(l+1)=brr(j) brr(j)=b crr(l+1)=crr(j) crr(j)=c jstack=jstack+2 if(jstack.gt.nstack)pause 'SORT2: NSTACK too small in sort2' if(ir-i+1.ge.j-l)then istack(jstack)=ir istack(jstack-1)=i ir=j-1 else istack(jstack)=j-1 istack(jstack-1)=l l=i endif endif goto 1 end C-----------------------------------------------------------------------