@@ -224,6 +224,38 @@ PURE REAL(RK) FUNCTION FREQUENCY_(WINDOW_TOTAL,WINDOW,SIGNAL)
224224 ! RESULT
225225 FREQUENCY_= (REAL (FST,RK)- 2._RK+2._RK * (FREQUENCY_-1._RK )/ REAL (NUM,RK))/ REAL (NUM,RK)
226226 END FUNCTION FREQUENCY_
227+
228+ PURE SUBROUTINE DECOMPOSITION_ (AVE ,WIN ,SIG ,ITR ,FRE ,COS_AMP ,SIN_AMP )
229+ ! $ACC ROUTINE SEQ
230+ REAL (RK),INTENT (IN ) :: AVE
231+ REAL (RK),DIMENSION (NUM),INTENT (IN ) :: WIN
232+ REAL (RK),DIMENSION (2 * NUM),INTENT (IN ) :: SIG
233+ INTEGER ,INTENT (IN ) :: ITR
234+ REAL (RK),DIMENSION (ITR),INTENT (OUT ) :: FRE
235+ REAL (RK),DIMENSION (ITR),INTENT (OUT ) :: COS_AMP
236+ REAL (RK),DIMENSION (ITR),INTENT (OUT ) :: SIN_AMP
237+ INTEGER :: I
238+ REAL (RK),DIMENSION (2 * NUM) :: ORB
239+ REAL (RK),DIMENSION (NUM) :: RAN
240+ REAL (RK),DIMENSION (2 * NUM) :: DEL
241+ ORB= SIG- AVE
242+ RAN= REAL ([(I,I= 1 ,NUM,1 )],RK)
243+ DO I= 1 ,ITR,1
244+ FRE(I)= FREQUENCY_(AVE,WIN,ORB)
245+ DEL(1 :2 * NUM:2 )=+ COS (TWO_PI* FRE(I)* RAN)
246+ DEL(2 :2 * NUM:2 )=- SIN (TWO_PI* FRE(I)* RAN)
247+ 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)
248+ 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)
249+ COS_AMP(I)= COS_AMP(I)/ REAL (NUM,RK)
250+ SIN_AMP(I)= SIN_AMP(I)/ REAL (NUM,RK)
251+ 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 )
252+ 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 )
253+ END DO
254+ IF (FLA== 0 )THEN
255+ COS_AMP = 2.0_RK * COS_AMP
256+ SIN_AMP = 2.0_RK * SIN_AMP
257+ END IF
258+ END SUBROUTINE DECOMPOSITION_
227259
228260END MODULE SIGNAL
229261
@@ -318,3 +350,91 @@ END MODULE SIGNAL
318350! !$ACC END DATA
319351! WRITE(*,*) OUT(1)-FRE
320352! END PROGRAM EXAMPLE
353+
354+ ! ! EXAMPLE (NUM=2**12,FLA=0)
355+ ! ! Real signal decomposition
356+ ! ! gfortran -o example -O3 -ffast-math -march=native signal.f90 example.f90
357+ ! ! example.f90:
358+ ! PROGRAM EXAMPLE
359+ ! USE SIGNAL
360+ ! IMPLICIT NONE
361+ ! REAL(RK),PARAMETER :: F1=1._RK*0.123456789_RK
362+ ! REAL(RK),PARAMETER :: F2=2._RK*0.123456789_RK
363+ ! REAL(RK),PARAMETER :: F3=3._RK*0.123456789_RK
364+ ! REAL(RK),PARAMETER :: F4=4._RK*0.123456789_RK
365+ ! REAL(RK),PARAMETER :: A1=1.0_RK
366+ ! REAL(RK),PARAMETER :: B1=1.E-1_RK
367+ ! REAL(RK),PARAMETER :: A2=1.E-1_RK
368+ ! REAL(RK),PARAMETER :: B2=5.E-2_RK
369+ ! REAL(RK),PARAMETER :: A3=1.E-3_RK
370+ ! REAL(RK),PARAMETER :: B3=0.0_RK
371+ ! REAL(RK),PARAMETER :: A4=0.0_RK
372+ ! REAL(RK),PARAMETER :: B4=1.E-4_RK
373+ ! REAL(RK),DIMENSION(2*NUM) :: SIG
374+ ! REAL(RK),DIMENSION(NUM) :: WIN
375+ ! REAL(RK) :: AVE
376+ ! INTEGER :: I
377+ ! REAL(RK),DIMENSION(4) :: FRE, COS_AMP, SIN_AMP
378+ ! SIG(1:2*NUM:2)=REAL([(I,I=1,NUM,1)],RK)
379+ ! 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)) + &
380+ ! A2*COS(TWO_PI*F2*SIG(1:2*NUM:2))+B2*SIN(TWO_PI*F2*SIG(1:2*NUM:2)) + &
381+ ! A3*COS(TWO_PI*F3*SIG(1:2*NUM:2))+B3*SIN(TWO_PI*F3*SIG(1:2*NUM:2)) + &
382+ ! A4*COS(TWO_PI*F4*SIG(1:2*NUM:2))+B4*SIN(TWO_PI*F4*SIG(1:2*NUM:2))
383+ ! SIG(2:2*NUM:2)=0.0_RK
384+ ! CALL WINDOW_(4,WIN)
385+ ! AVE = SUM(WIN)
386+ ! FRE = 0.0_RK
387+ ! COS_AMP=0.0_RK
388+ ! SIN_AMP=0.0_RK
389+ ! CALL DECOMPOSITION_(AVE,WIN,SIG,4,FRE,COS_AMP,SIN_AMP)
390+ ! DO I=1,4,1
391+ ! WRITE(*,*) FRE(I),COS_AMP(I),SIN_AMP(I)
392+ ! END DO
393+ ! END PROGRAM EXAMPLE
394+
395+
396+ ! ! EXAMPLE (NUM=2**12,FLA=1)
397+ ! ! Complex signal decomposition
398+ ! ! gfortran -o example -O3 -ffast-math -march=native signal.f90 example.f90
399+ ! ! example.f90:
400+ ! PROGRAM EXAMPLE
401+ ! USE SIGNAL
402+ ! IMPLICIT NONE
403+ ! REAL(RK),PARAMETER :: F1=1._RK*0.123456789_RK
404+ ! REAL(RK),PARAMETER :: F2=2._RK*0.123456789_RK
405+ ! REAL(RK),PARAMETER :: F3=3._RK*0.123456789_RK
406+ ! REAL(RK),PARAMETER :: F4=4._RK*0.123456789_RK
407+ ! REAL(RK),PARAMETER :: A1=1.0_RK
408+ ! REAL(RK),PARAMETER :: B1=1.E-1_RK
409+ ! REAL(RK),PARAMETER :: A2=1.E-1_RK
410+ ! REAL(RK),PARAMETER :: B2=5.E-2_RK
411+ ! REAL(RK),PARAMETER :: A3=1.E-3_RK
412+ ! REAL(RK),PARAMETER :: B3=0.0_RK
413+ ! REAL(RK),PARAMETER :: A4=0.0_RK
414+ ! REAL(RK),PARAMETER :: B4=1.E-4_RK
415+ ! REAL(RK),DIMENSION(2*NUM) :: SIG
416+ ! REAL(RK),DIMENSION(NUM) :: WIN
417+ ! REAL(RK) :: AVE
418+ ! INTEGER :: I
419+ ! REAL(RK),DIMENSION(4) :: FRE, COS_AMP, SIN_AMP
420+ ! SIG(1:2*NUM:2)=REAL([(I,I=1,NUM,1)],RK)
421+ ! 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)) + &
422+ ! A2*COS(TWO_PI*F2*SIG(1:2*NUM:2))+B2*SIN(TWO_PI*F2*SIG(1:2*NUM:2)) + &
423+ ! A3*COS(TWO_PI*F3*SIG(1:2*NUM:2))+B3*SIN(TWO_PI*F3*SIG(1:2*NUM:2)) + &
424+ ! A4*COS(TWO_PI*F4*SIG(1:2*NUM:2))+B4*SIN(TWO_PI*F4*SIG(1:2*NUM:2))
425+ !
426+ ! SIG(2:2*NUM:2)=REAL([(I,I=1,NUM,1)],RK)
427+ ! 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)) + &
428+ ! B2*COS(TWO_PI*F2*SIG(2:2*NUM:2))-A2*SIN(TWO_PI*F2*SIG(2:2*NUM:2)) + &
429+ ! B3*COS(TWO_PI*F3*SIG(2:2*NUM:2))-A3*SIN(TWO_PI*F3*SIG(2:2*NUM:2)) + &
430+ ! B4*COS(TWO_PI*F4*SIG(2:2*NUM:2))-A4*SIN(TWO_PI*F4*SIG(2:2*NUM:2))
431+ ! CALL WINDOW_(4,WIN)
432+ ! AVE = SUM(WIN)
433+ ! FRE = 0.0_RK
434+ ! COS_AMP=0.0_RK
435+ ! SIN_AMP=0.0_RK
436+ ! CALL DECOMPOSITION_(AVE,WIN,SIG,4,FRE,COS_AMP,SIN_AMP)
437+ ! DO I=1,4,1
438+ ! WRITE(*,*) FRE(I),COS_AMP(I),SIN_AMP(I)
439+ ! END DO
440+ ! END PROGRAM EXAMPLE
0 commit comments