1. 程式人生 > >fortran實現正餘弦變換

fortran實現正餘弦變換

    本來是想正餘弦變換都做的,根據前人研究:正弦變換的結果比餘弦更精確。所以本程式碼只做了正弦變化,如果讀者想實現餘弦變換,只需要把對應的正弦變換系數改為餘弦變換系數即可。在末尾本文會給出250點正餘弦變換系數。
    程式碼如下:
!.. 本程式碼將實現正餘弦變換
!.. 測試函式
!.. ∫exp(-w)sin(wt)dw
!.. 正弦變換離散表示式
!.. f(t)=∑{K(w)*wsl} * coef
!.. coef = 1/t*sqrt(pi/2)
Module testSineTransform
  implicit none
  integer, parameter :: nt = 46
  real(kind=8), parameter :: pi = acos(-1.d0)
  real(kind=8) :: t(nt), t1 = 5, t2 = 50 !.. 單位:s
  real(kind=8) :: wsl(-149:100), w, delta
contains
  subroutine cal_SineTransform
    implicit none
    integer :: i, j, fileid
    real(kind=8) :: s1(nt), s2(nt) !.. s1為數值解, s2為解析解

    do i = 1, nt
        t(i) = t1 + dble(i-1) * ( t2 - t1 ) / dble( nt - 1)
        s2(i) = func2( t(i) )!.. 解析表示式
    end do

    open( newunit = fileid, file = "sin250.txt" )
    do i = -149, 100
        read( fileid,* ) wsl(i)
    end do
    close( fileid )

    !.. 計算正弦變換
    delta = log(10.d0) / 20.d0 !.. 理論推導而來,無需更改
    s1 = 0.d0

    open( newunit = fileid, file = "result.dat" )
    do i = 1, nt
        do j = -149, 100
            w = exp( dble(j)*delta ) / t(i) !.. 理論推導而來,無需更改
            s1(i) = s1(i) + func1( w ) * wsl(j)
        end do
    s1(i) = s1(i) / t(i) * sqrt( pi/2.d0 ) !.. 1/t*sqrt(pi/2)為正弦變換出來的係數
    write( fileid, '(*(g0,3x))' ) t(i), s1(i), s2(i), abs( s1(i)-s2(i) )/s2(i)*100
    end do
    close( fileid )
  end subroutine cal_SineTransform

  real(kind=8) function func1(w) !.. 核函式
    implicit none
    real(kind=8) :: w

    func1 = exp(-1.d0*w)

  end function func1

  real(kind=8) function func2(t) !.. 解析解
    implicit none
    real(kind=8) :: t

    func2 = t / ( t*t + 1.d0 )

  end function func2
End module testSineTransform

Program main
  use testSineTransform
  implicit none
  call cal_SineTransform
End program main

下面給出250點正弦變換系數
正弦變換系數:

 1.156447054648637e-016
 1.455880584491685e-016
 1.832845064354326e-016
 2.307415227397174e-016
 2.904863665331116e-016
 3.657006686082886e-016
 4.603898648210914e-016
 5.795965001557580e-016
 7.296687626330448e-016
 9.185985474711453e-016
 1.156447054648634e-015
 1.455880584491681e-015
 1.832845064354321e-015
 2.307415227397166e-015
 2.904863665331103e-015
 3.657006686082864e-015
 4.603898648210901e-015
 5.795965001557501e-015
 7.296687626330329e-015
 9.185985474711278e-015
 1.156447054648612e-014
 1.455880584491647e-014
 1.832845064354265e-014
 2.307415227397080e-014
 2.904863665330965e-014
 3.657006686082646e-014
 4.603898648210535e-014
 5.795965001556953e-014
 7.296687626329462e-014
 9.185985474709900e-014
 1.156447054648389e-013
 1.455880584491301e-013
 1.832845064353717e-013
 2.307415227396208e-013
 2.904863665329599e-013
 3.657006686080477e-013
 4.603898648207092e-013
 5.795965001551463e-013
 7.296687626320796e-013
 9.185985474696154e-013
 1.156447054646209e-012
 1.455880584487839e-012
 1.832845064348231e-012
 2.307415227387515e-012
 2.904863665315808e-012
 3.657006686058640e-012
 4.603898648172482e-012
 5.795965001496609e-012
 7.296687626233827e-012
 9.185985474558366e-012
 1.156447054624371e-011
 1.455880584453228e-011
 1.832845064293375e-011
 2.307415227300576e-011
 2.904863665178019e-011
 3.657006685840241e-011
 4.603898647826370e-011
 5.795965000948059e-011
 7.296687625364435e-011
 9.185985473180428e-011
 1.156447054405990e-010
 1.455880584107116e-010
 1.832845063744824e-010
 2.307415226431180e-010
 2.904863663800132e-010
 3.657006683656435e-010
 4.603898644365244e-010
 5.795964995462562e-010
 7.296687616670500e-010
 9.185985459401473e-010
 1.156447052222165e-009
 1.455880580645994e-009
 1.832845058259313e-009
 2.307415217737235e-009
 2.904863650021138e-009
 3.657006661818234e-009
 4.603898609754006e-009
 5.795964940607466e-009
 7.296687529731007e-009
 9.185985321611760e-009
 1.156447030383951e-008
 1.455880546034779e-008
 1.832845003404190e-008
 2.307415130797806e-008
 2.904863512231297e-008
 3.657006443436261e-008
 4.603898263641547e-008
 5.795964392056795e-008
 7.296686660335758e-008
 9.185983943714872e-008
 1.156446812001731e-007
 1.455880199922772e-007
 1.832844454852834e-007
 2.307414261403871e-007
 2.904862134332549e-007
 3.657004259617825e-007
 4.603894802516352e-007
 5.795958906554214e-007
 7.296677966382738e-007
 9.185970164759857e-007
 1.156444628179741e-006
 1.455876738807346e-006
 1.832838969341425e-006
 2.307405567480775e-006
 2.904848355357452e-006
 3.656982421491810e-006
 4.603860191324607e-006
 5.795904051742752e-006
 7.296591027122000e-006
 9.185832376013115e-006
 1.156422790075531e-005
 1.455842127958865e-005
 1.832784114711636e-005
 2.307318629427335e-005
 2.904710567601305e-005
 3.656764044807223e-005
 4.603514087536846e-005
 5.795355521527330e-005
 7.295721667410295e-005
 9.184454558911375e-005
 1.156204422278983e-004
 1.455496047235377e-004
 1.832235621566062e-004
 2.306449358701818e-004
 2.903332902825535e-004
 3.654580713105382e-004
 4.600053900083832e-004
 5.789871945981126e-004
 7.287031465075219e-004
 9.170683245122596e-004
 1.154022098295481e-003
 1.452037965023791e-003
 1.826756086858394e-003
 2.297767484717661e-003
 2.889577754406488e-003
 3.632790471922652e-003
 4.565537633632412e-003
 5.735207473181395e-003
 7.200470157233521e-003
 9.033651050110031e-003
 1.132334380467910e-002
 1.417727733464931e-002
 1.772499346509631e-002
 2.212023128150156e-002
 2.754164273887034e-002
 3.419148067351661e-002
 4.228846080133485e-002
 5.205416439533162e-002
 6.368324102630288e-002
 7.729769668970302e-002
 9.286191525747088e-002
 1.100669688846222e-001
 1.281325720732841e-001
 1.455759302997095e-001
 1.598450602566018e-001
 1.675994998118554e-001
 1.619683745893242e-001
 1.351613440916820e-001
 8.220710085812771e-002
-5.378512256133189e-003
-1.293051042068928e-001
-2.723112991886841e-001
-3.978468461558793e-001
-4.394993669496010e-001
-3.176425577428793e-001
 1.665756359292827e-002
 4.630044483056220e-001
 7.298213861736705e-001
 3.957559489684660e-001
-4.902404478827365e-001
-1.020512348109581e+000
 5.251407980261734e-002
 1.253864451347815e+000
-1.317577223583377e-001
-1.549527153893942e+000
 1.822344193602740e+000
-1.058204213789976e+000
 2.629725097214586e-001
 1.625578365785221e-001
-2.871269632839427e-001
 2.688993168200287e-001
-2.085718537012279e-001
 1.489333044554680e-001
-1.023510671993418e-001
 6.920616476794303e-002
-4.650126535306088e-002
 3.114481903045837e-002
-2.081379273802720e-002
 1.390327480927853e-002
-9.298886523349949e-003
 6.224588805575879e-003
-4.160927227996026e-003
 2.775386859464604e-003
-1.851896921274413e-003
 1.239790988028041e-003
-8.315445781273338e-004
 5.558587783122442e-004
-3.696533331522879e-004
 2.460602463017792e-004
-1.651618293595525e-004
 1.114024919622099e-004
-7.389805836530149e-005
 4.938931786488272e-005
-3.300905019751596e-005
 2.206139792679832e-005
-1.474460111811994e-005
 9.854464474452981e-006
-6.586171392444854e-006
 4.401827590217889e-006
-2.941934696111318e-006
 1.966224160030928e-006
-1.314113957933395e-006
 8.782800707768828e-007
-5.869931432254631e-007
 3.923132969292148e-007
-2.622002057839472e-007
 1.752399127209358e-007
-1.171205297822945e-007
 7.827679370240814e-008
-5.231581895777865e-008
 3.496495939305279e-008
-2.336861793837973e-008
 1.561827366110047e-008
-1.043837820431878e-008
 6.976426582137068e-009
-4.662652272535390e-009
 3.116255286086561e-009
-2.082730266047119e-009
 1.391980105248554e-009
-9.303214367194074e-010
 6.217746736151215e-010
-4.155593212088843e-010
 2.777365447188800e-010
-1.856235303493902e-010
 1.240603574666303e-010
-8.291498532420789e-011
 5.541572611672626e-011
-3.703676348776440e-011
 2.475329560347511e-011
-1.654371455635934e-011
 1.105689099773430e-011
-7.389805845555023e-012
 4.938931788889779e-012
-3.300905020390523e-012
 2.206139792849932e-012
-1.474460111857239e-012
 9.854464474573356e-013
-6.586171392476662e-013
 4.401827590226633e-013
-2.941934696113489e-013


餘弦變換系數:

250點餘弦變換濾波係數(王華軍, 2004)
 n=-149~-100	  n=-99~-50	  n=-49~0	       n=1~50	      n=51~100
3.25931064E-09	1.03068452E-06	3.25929012E-04	4.47184906E-02	-1.61220740E-04
3.65700668E-09	1.15644705E-06	3.65697770E-04	3.54831893E-02	1.08092999E-04
4.10322899E-09	1.29755494E-06	4.10318805E-04	2.04513710E-02	-7.22432985E-05
4.60389865E-09	1.45588058E-06	4.60384082E-04	-2.05568659E-03	4.82833692E-05
5.16565924E-09	1.63352488E-06	5.16557757E-04	-3.36444840E-02	-3.22698961E-05
5.79596500E-09	1.83284506E-06	5.79584963E-04	-7.63113734E-02	2.15673888E-05
6.50317969E-09	2.05648598E-06	6.50301672E-04	-1.26643432E-01	-1.44144331E-05
7.29668762E-09	2.30741523E-06	7.29645743E-04	-1.86274938E-01	9.63379869E-06
8.18701817E-09	2.58896246E-06	8.18669301E-04	-2.46135314E-01	-6.43869076E-06
9.18598547E-09	2.90486366E-06	9.18552617E-04	-2.90298367E-01	4.30325981E-06
1.03068452E-08	3.25931064E-06	1.03061964E-03	-2.98807588E-01	-2.87605752E-06
1.15644705E-08	3.65700668E-06	1.15635541E-03	-2.44850614E-01	1.92219555E-06
1.29755494E-08	4.10322898E-06	1.29742549E-03	-9.95001792E-02	-1.28468770E-06
1.45588058E-08	4.60389864E-06	1.45569774E-03	1.35603930E-01	8.58613212E-07
1.63352488E-08	5.16565923E-06	1.63326661E-03	4.09054093E-01	-5.73848918E-07
1.83284506E-08	5.79596499E-06	1.83248024E-03	5.77439580E-01	3.83528434E-07
2.05648598E-08	6.50317967E-06	2.05597066E-03	4.57925568E-01	-2.56328896E-07
2.30741523E-08	7.29668760E-06	2.30668732E-03	-6.98345034E-02	1.71315859E-07
2.58896247E-08	8.18701814E-06	2.58793429E-03	-7.05728308E-01	-1.14497912E-07
2.90486366E-08	9.18598542E-06	2.90341135E-03	-7.84315028E-01	7.65239831E-08
3.25931064E-08	1.03068452E-05	3.25725924E-03	2.52561836E-01	-5.11443386E-08
3.65700668E-08	1.15644704E-05	3.65410908E-03	1.11765483E+00	3.41820077E-08
4.10322899E-08	1.29755492E-05	4.09913616E-03	7.96337157E-02	-2.28453369E-08
4.60389865E-08	1.45588057E-05	4.59811762E-03	-1.54929662E+00	1.52685420E-08
5.16565924E-08	1.63352486E-05	5.15749377E-03	9.69765096E-01	-1.02046372E-08
5.79596500E-08	1.83284503E-05	5.78443175E-03	5.22594045E-01	6.82020717E-09
6.50317969E-08	2.05648593E-05	6.48688994E-03	-1.38177462E+00	-4.55824398E-09
7.29668762E-08	2.30741515E-05	7.27368023E-03	1.42319844E+00	3.04647464E-09
8.18701817E-08	2.58896236E-05	8.15452380E-03	-1.10540419E+00	-2.03609279E-09
9.18598547E-08	2.90486352E-05	9.14009379E-03	7.66503550E-01	1.36081024E-09
1.03068452E-07	3.25931043E-05	1.02420355E-02	-5.12775503E-01	-9.09489244E-10
1.15644705E-07	3.65700639E-05	1.14729492E-02	3.41609163E-01	6.07851604E-10
1.29755494E-07	4.10322858E-05	1.28463164E-02	-2.28175614E-01	-4.06253922E-10
1.45588058E-07	4.60389807E-05	1.43763378E-02	1.52559552E-01	2.71517338E-10
1.63352488E-07	5.16565843E-05	1.60776460E-02	-1.01971829E-01	-1.81466961E-10
1.83284506E-07	5.79596384E-05	1.79648247E-02	6.81663603E-02	1.21282340E-10
2.05648598E-07	6.50317806E-05	2.00516669E-02	-4.55882998E-02	-8.10583142E-11
2.30741523E-07	7.29668532E-05	2.23500306E-02	3.04810683E-02	5.41748312E-11
2.58896247E-07	8.18701492E-05	2.48681715E-02	-2.03615795E-02	-3.62074189E-11
2.90486366E-07	9.18598088E-05	2.76082657E-02	1.35977270E-02	2.41990082E-11
3.25931064E-07	1.03068387E-04	3.05629294E-02	-9.08987662E-03	-1.61732600E-11
3.65700668E-07	1.15644614E-04	3.37101508E-02	6.08272853E-03	1.08092999E-11
4.10322899E-07	1.29755364E-04	3.70064047E-02	-4.06761271E-03	-7.22432985E-12
4.60389865E-07	1.45587876E-04	4.03767590E-02	2.71487245E-03	4.82833692E-12
5.16565924E-07	1.63352230E-04	4.37019801E-02	-1.81128418E-03	-3.22698961E-12
5.79596500E-07	1.83284142E-04	4.68002258E-02	1.21137650E-03	2.15673888E-12
6.50317969E-07	2.05648083E-04	4.94047004E-02	-8.12064069E-04	-1.44144331E-12
7.29668762E-07	2.30740795E-04	5.11324187E-02	5.43426304E-04	9.63379869E-13
8.18701817E-07	2.58895218E-04	5.14507096E-02	-3.61946442E-04	-6.43869076E-13
9.18598547E-07	2.90484914E-04	4.96316608E-02	2.40820925E-04	4.30325981E-13