Skip to content

Commit da70ae2

Browse files
author
Ivan Morozov
committed
FORTRAN functions (+decomposition)
1 parent df9d72e commit da70ae2

File tree

1 file changed

+124
-2
lines changed

1 file changed

+124
-2
lines changed

‎signal.f90‎

Lines changed: 124 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -19,8 +19,8 @@ MODULE SIGNAL
1919

2020
PRIVATE
2121

22-
INTEGER,PARAMETER :: RK=REAL64 ! REAL KIND
23-
INTEGER,PARAMETER :: NUM=1024 ! SIGNAL LENGTH (POWER OF TWO)
22+
INTEGER,PARAMETER :: RK=REAL64 ! REAL KIND
23+
INTEGER,PARAMETER :: NUM=4096 ! SIGNAL LENGTH (POWER OF TWO)
2424
INTEGER,PARAMETER :: FLA=1 ! COMPLEX SIGNAL FLAG (0/1 FOR REAL/COMPLEX INPUT)
2525
REAL(RK),PARAMETER :: PI=3.141592653589793238460_RK ! PI
2626
REAL(RK),PARAMETER :: TWO_PI=2.0_RK*PI ! 2*PI
@@ -38,6 +38,7 @@ MODULE SIGNAL
3838
PUBLIC :: TWO_PI
3939
PUBLIC :: WINDOW_
4040
PUBLIC :: FREQUENCY_
41+
PUBLIC :: DECOMPOSITION_
4142

4243
PUBLIC :: BRENT_MAX
4344
PUBLIC :: BRENT_TOL
@@ -247,6 +248,38 @@ PURE REAL(RK) FUNCTION FREQUENCY_(WINDOW_TOTAL,WINDOW,SIGNAL)
247248
FREQUENCY_=(REAL(FST,RK)-2._RK+2._RK*(FREQUENCY_-1._RK)/REAL(NUM,RK))/REAL(NUM,RK)
248249
END FUNCTION FREQUENCY_
249250

251+
PURE SUBROUTINE DECOMPOSITION_(AVE,WIN,SIG,ITR,FRE,COS_AMP,SIN_AMP)
252+
!$ACC ROUTINE SEQ
253+
REAL(RK),INTENT(IN) :: AVE
254+
REAL(RK),DIMENSION(NUM),INTENT(IN) :: WIN
255+
REAL(RK),DIMENSION(2*NUM),INTENT(IN) :: SIG
256+
INTEGER,INTENT(IN) :: ITR
257+
REAL(RK),DIMENSION(ITR),INTENT(OUT) :: FRE
258+
REAL(RK),DIMENSION(ITR),INTENT(OUT) :: COS_AMP
259+
REAL(RK),DIMENSION(ITR),INTENT(OUT) :: SIN_AMP
260+
INTEGER :: I
261+
REAL(RK),DIMENSION(2*NUM) :: ORB
262+
REAL(RK),DIMENSION(NUM) :: RAN
263+
REAL(RK),DIMENSION(2*NUM) :: DEL
264+
ORB=SIG-AVE
265+
RAN=REAL([(I,I=1,NUM,1)],RK)
266+
DO I=1,ITR,1
267+
FRE(I)=FREQUENCY_(AVE,WIN,ORB)
268+
DEL(1:2*NUM:2)=+COS(TWO_PI*FRE(I)*RAN)
269+
DEL(2:2*NUM:2)=-SIN(TWO_PI*FRE(I)*RAN)
270+
COS_AMP(I)=SUM(SIG(1:2*NUM:2)*DEL(1:2*NUM:2)*WIN)+SUM(SIG(2:2*NUM:2)*DEL(2:2*NUM:2)*WIN)
271+
SIN_AMP(I)=SUM(SIG(2:2*NUM:2)*DEL(1:2*NUM:2)*WIN)-SUM(SIG(1:2*NUM:2)*DEL(2:2*NUM:2)*WIN)
272+
COS_AMP(I)=COS_AMP(I)/REAL(NUM,RK)
273+
SIN_AMP(I)=SIN_AMP(I)/REAL(NUM,RK)
274+
ORB(1:2*NUM:2)=ORB(1:2*NUM:2)-COS_AMP(I)*DEL(1:2*NUM:2)+SIN_AMP(I)*DEL(2:2*NUM:2)
275+
ORB(2:2*NUM:2)=ORB(2:2*NUM:2)-COS_AMP(I)*DEL(2:2*NUM:2)-SIN_AMP(I)*DEL(1:2*NUM:2)
276+
END DO
277+
IF(FLA==0)THEN
278+
COS_AMP = 2.0_RK*COS_AMP
279+
SIN_AMP = 2.0_RK*SIN_AMP
280+
END IF
281+
END SUBROUTINE DECOMPOSITION_
282+
250283
! POLINT (NRF90)
251284
PURE SUBROUTINE POLINT_(XA,YA,X,Y,DY)
252285
REAL(RK),DIMENSION(:),INTENT(IN) :: XA
@@ -559,3 +592,92 @@ END MODULE SIGNAL
559592
! !$ACC END DATA
560593
! WRITE(*,*) OUT(1)-FRE
561594
! END PROGRAM EXAMPLE
595+
596+
! ! EXAMPLE (NUM=2**12,FLA=0)
597+
! ! Real signal decomposition
598+
! ! gfortran -o example -O3 -ffast-math -march=native signal.f90 example.f90
599+
! ! example.f90:
600+
! PROGRAM EXAMPLE
601+
! USE SIGNAL
602+
! IMPLICIT NONE
603+
! REAL(RK),PARAMETER :: F1=1._RK*0.123456789_RK
604+
! REAL(RK),PARAMETER :: F2=2._RK*0.123456789_RK
605+
! REAL(RK),PARAMETER :: F3=3._RK*0.123456789_RK
606+
! REAL(RK),PARAMETER :: F4=4._RK*0.123456789_RK
607+
! REAL(RK),PARAMETER :: A1=1.0_RK
608+
! REAL(RK),PARAMETER :: B1=1.E-1_RK
609+
! REAL(RK),PARAMETER :: A2=1.E-1_RK
610+
! REAL(RK),PARAMETER :: B2=5.E-2_RK
611+
! REAL(RK),PARAMETER :: A3=1.E-3_RK
612+
! REAL(RK),PARAMETER :: B3=0.0_RK
613+
! REAL(RK),PARAMETER :: A4=0.0_RK
614+
! REAL(RK),PARAMETER :: B4=1.E-4_RK
615+
! REAL(RK),DIMENSION(2*NUM) :: SIG
616+
! REAL(RK),DIMENSION(NUM) :: WIN
617+
! REAL(RK) :: AVE
618+
! INTEGER :: I
619+
! REAL(RK),DIMENSION(4) :: FRE, COS_AMP, SIN_AMP
620+
! SIG(1:2*NUM:2)=REAL([(I,I=1,NUM,1)],RK)
621+
! SIG(1:2*NUM:2)=A1*COS(TWO_PI*F1*SIG(1:2*NUM:2))+B1*SIN(TWO_PI*F1*SIG(1:2*NUM:2)) + &
622+
! A2*COS(TWO_PI*F2*SIG(1:2*NUM:2))+B2*SIN(TWO_PI*F2*SIG(1:2*NUM:2)) + &
623+
! A3*COS(TWO_PI*F3*SIG(1:2*NUM:2))+B3*SIN(TWO_PI*F3*SIG(1:2*NUM:2)) + &
624+
! A4*COS(TWO_PI*F4*SIG(1:2*NUM:2))+B4*SIN(TWO_PI*F4*SIG(1:2*NUM:2))
625+
! SIG(2:2*NUM:2)=0.0_RK
626+
! CALL WINDOW_(4,WIN)
627+
! AVE = SUM(WIN)
628+
! FRE = 0.0_RK
629+
! COS_AMP=0.0_RK
630+
! SIN_AMP=0.0_RK
631+
! CALL DECOMPOSITION_(AVE,WIN,SIG,4,FRE,COS_AMP,SIN_AMP)
632+
! DO I=1,4,1
633+
! WRITE(*,*) FRE(I),COS_AMP(I),SIN_AMP(I)
634+
! END DO
635+
! END PROGRAM EXAMPLE
636+
637+
638+
! ! EXAMPLE (NUM=2**12,FLA=1)
639+
! ! Complex signal decomposition
640+
! ! gfortran -o example -O3 -ffast-math -march=native signal.f90 example.f90
641+
! ! example.f90:
642+
! PROGRAM EXAMPLE
643+
! USE SIGNAL
644+
! IMPLICIT NONE
645+
! REAL(RK),PARAMETER :: F1=1._RK*0.123456789_RK
646+
! REAL(RK),PARAMETER :: F2=2._RK*0.123456789_RK
647+
! REAL(RK),PARAMETER :: F3=3._RK*0.123456789_RK
648+
! REAL(RK),PARAMETER :: F4=4._RK*0.123456789_RK
649+
! REAL(RK),PARAMETER :: A1=1.0_RK
650+
! REAL(RK),PARAMETER :: B1=1.E-1_RK
651+
! REAL(RK),PARAMETER :: A2=1.E-1_RK
652+
! REAL(RK),PARAMETER :: B2=5.E-2_RK
653+
! REAL(RK),PARAMETER :: A3=1.E-3_RK
654+
! REAL(RK),PARAMETER :: B3=0.0_RK
655+
! REAL(RK),PARAMETER :: A4=0.0_RK
656+
! REAL(RK),PARAMETER :: B4=1.E-4_RK
657+
! REAL(RK),DIMENSION(2*NUM) :: SIG
658+
! REAL(RK),DIMENSION(NUM) :: WIN
659+
! REAL(RK) :: AVE
660+
! INTEGER :: I
661+
! REAL(RK),DIMENSION(4) :: FRE, COS_AMP, SIN_AMP
662+
! SIG(1:2*NUM:2)=REAL([(I,I=1,NUM,1)],RK)
663+
! SIG(1:2*NUM:2)=A1*COS(TWO_PI*F1*SIG(1:2*NUM:2))+B1*SIN(TWO_PI*F1*SIG(1:2*NUM:2)) + &
664+
! A2*COS(TWO_PI*F2*SIG(1:2*NUM:2))+B2*SIN(TWO_PI*F2*SIG(1:2*NUM:2)) + &
665+
! A3*COS(TWO_PI*F3*SIG(1:2*NUM:2))+B3*SIN(TWO_PI*F3*SIG(1:2*NUM:2)) + &
666+
! A4*COS(TWO_PI*F4*SIG(1:2*NUM:2))+B4*SIN(TWO_PI*F4*SIG(1:2*NUM:2))
667+
!
668+
! SIG(2:2*NUM:2)=REAL([(I,I=1,NUM,1)],RK)
669+
! SIG(2:2*NUM:2)=B1*COS(TWO_PI*F1*SIG(2:2*NUM:2))-A1*SIN(TWO_PI*F1*SIG(2:2*NUM:2)) + &
670+
! B2*COS(TWO_PI*F2*SIG(2:2*NUM:2))-A2*SIN(TWO_PI*F2*SIG(2:2*NUM:2)) + &
671+
! B3*COS(TWO_PI*F3*SIG(2:2*NUM:2))-A3*SIN(TWO_PI*F3*SIG(2:2*NUM:2)) + &
672+
! B4*COS(TWO_PI*F4*SIG(2:2*NUM:2))-A4*SIN(TWO_PI*F4*SIG(2:2*NUM:2))
673+
! CALL WINDOW_(4,WIN)
674+
! AVE = SUM(WIN)
675+
! FRE = 0.0_RK
676+
! COS_AMP=0.0_RK
677+
! SIN_AMP=0.0_RK
678+
! CALL DECOMPOSITION_(AVE,WIN,SIG,4,FRE,COS_AMP,SIN_AMP)
679+
! DO I=1,4,1
680+
! WRITE(*,*) FRE(I),COS_AMP(I),SIN_AMP(I)
681+
! END DO
682+
! END PROGRAM EXAMPLE
683+

0 commit comments

Comments
 (0)