none
Obtainer of THD-F, THD-R, THD+N, Vrms RRS feed

  • Question


  • Function GetNumF1(ByVal f1 As Double, ByVal Numb As Long, ByVal fs As Double) As Long

     GetNumF1 = 1 + Int(f1 * Numb / fs)

    End Function

    Function GetTHD_F(ByRef xn() As Double, ByVal Numb As Long, ByVal fs As Double, ByVal f1 As Double) As Double


    Dim numbf1 As Long
    Dim ModXk(1 To NMAX) As Double
    Dim argXk(1 To NMAX) As Double
    Dim Urms1 As Double
    Dim Uhharm   As Double
    Dim i As Long


    numbf1 = Int(f1 * Numb / fs) ' 2  , get associated vector of frequency, ORIGIN=0
    Call GetSpectrum(xn, Numb, ModXk, argXk, 0)
    Urms1 = ModXk(numbf1 + 1) / ((2) ^ (1 / 2)) ' / 1.41

    i = numbf1 * 2
    Uhharm = 0
    Do While (i <= Numb + 1)
     If (i + 1 <= Numb) Then
     Uhharm = Uhharm + ((ModXk(i + 1) ^ 2) / 2)    ' 2*2 +1, 3*2+1
     End If
    i = i + numbf1
    Loop

    Uhharm = (Uhharm) ^ (1 / 2)
    GetTHD_F = Uhharm * 100 / (Urms1)

    End Function


    Function GetTHD_R(ByRef xn() As Double, ByVal Numb As Long, ByVal fs As Double, ByVal f1 As Double) As Double

    Dim numbf1 As Long
    Dim ModXk(1 To NMAX) As Double
    Dim argXk(1 To NMAX) As Double
    Dim Urmssigah As Double
    Dim Uhharm   As Double
    Dim i As Long


     
    numbf1 = Int(f1 * Numb / fs)    ' 2 , get associated vector of frequency, ORIGIN=0
     
    Call GetSpectrum(xn, Numb, ModXk, argXk, 0)

    Urmssigah = 0
    i = numbf1

    Do While (i <= Numb)

    If (i + 1 <= Numb) Then
    Urmssigah = Urmssigah + (ModXk(i + 1) ^ 2) / 2
    End If
    i = i + numbf1
    Loop

    Urmssigah = (Urmssigah) ^ (0.5)

    Uhharm = 0


    ' 2*2=4

    i = numbf1 * 2 ' number of  bin of 2-nd harmonic , for example, 4
    Do While (i <= Numb)
                            ' ModXk(4+1=5)
    If (i + 1 <= Numb) Then
    Uhharm = Uhharm + (ModXk(i + 1) ^ 2) / 2
    End If
    i = i + numbf1          ' 4+2=6,next harmonic= ModXk(6+1=7)
    Loop
    Uhharm = (Uhharm) ^ (0.5)

    GetTHD_R = Uhharm * 100 / Urmssigah


    End Function

     Function GetRMS(ByRef xn() As Double, ByVal Numb As Long)
    Dim i As Long
    Dim Urms As Double
    Urms = 0
    For i = 1 To Numb
     Urms = Urms + xn(i) ^ 2
    Next i
    Urms = Urms / Numb
    GetRMS = (Urms) ^ (0.5)


    End Function


    Function GetTHD_N(ByRef xn() As Double, ByVal Numb As Long, ByVal fs As Double, ByVal f1 As Double)

    Dim numbf1 As Long
    Dim ModXk(1 To NMAX) As Double
    Dim argXk(1 To NMAX) As Double
    Dim Urmssig As Double
    Dim Uhharmn   As Double
    Dim i As Long

     

    Urmssig = 0
    For i = 1 To Numb
     Urmssig = Urmssig + xn(i) ^ 2
    Next i
    Urmssig = Urmssig / Numb
    Urmssig = (Urmssig) ^ (0.5)
    'Cells(5, 1).Value = Urmssig

     
    numbf1 = 1 + Int(f1 * Numb / fs)
    'Cells(6, 1).Value = numbf1

    Call GetSpectrum(xn, Numb, ModXk, argXk, 0)

    Uhharmn = 0
    i = 1

    Do While (i <= Numb)
    Cells(i + 1, 7) = ModXk(i)
    If (i <> numbf1) Then


    Uhharmn = Uhharmn + 0.5 * (ModXk(i) ^ 2)

    End If

    i = i + 1

    Loop
    Uhharmn = (Uhharmn) ^ (0.5)
    'Cells(13, 1) = Uhharmn
    GetTHD_N = Uhharmn * 100 / Urmssig


    End Function


    Sub TestTHDN()
    Dim modX(1 To NMAX) As Double
    Dim argX(1 To NMAX) As Double
    Dim Xi(1 To NMAX)   As Double
    Dim Xq(1 To NMAX)   As Double
    Dim fvect(1 To NMAX) As Double
    Dim fs As Double
    Dim f1 As Double
    Dim N As Long
    Dim i As Long

    N = Cells(1, 1).Value
    fs = Cells(2, 1).Value
    f1 = Cells(3, 1).Value
    Call Getfassoc(fs, fvect, N)

    For i = 1 To N

     Cells(i + 1, 2).Value = fvect(i)

    modX(i) = Cells(i + 1, 3).Value
    argX(i) = Cells(i + 1, 4).Value

    Next i


      Call SynthSignalfromRFi(modX, argX, Xi, Xq, N, 1)
      Call printcomplex(Xi, Xq, N, 5, 6)
     
     
     Cells(15, 1).Value = GetTHD_N(Xi, N, fs, f1)
     Cells(16, 1).Value = GetTHD_R(Xi, N, fs, f1)
     Cells(17, 1).Value = GetTHD_F(Xi, N, fs, f1)


    End Sub

    Friday, November 27, 2015 8:45 PM

All replies

  • Source with getspectrum, synth signal

     Option Explicit
     Const M_PI As Double = 3.14159265358979
     
      Const NMAX As Long = 131072

     Sub Getfassoc(ByVal fs As Double, ByRef fvect() As Double, ByVal Numb As Long)
       Dim k As Long
        For k = 0 To Numb - 1
         fvect(k + 1) = fs * k / Numb
        Next k
      End Sub
     
     
      Sub Gettassoc(ByVal fs As Double, ByRef tvect() As Double, ByVal Numb As Long)
       Dim N As Long
        For N = 0 To Numb - 1
         tvect(N + 1) = N / fs
        Next N
      End Sub


      Function Arg(ByVal x_re As Double, ByVal x_im As Double) As Double
       'Dim x_re,x_im As Double
       'Dim arg As Double
       Dim message As String
      If (x_re > 0) Then
                 Arg = Atn(x_im / x_re)
                 End If
      If (x_re = 0) And (x_im > 0) Then
                 Arg = M_PI / 2
                 End If
      If (x_re = 0) And (x_im < 0) Then
                 Arg = -M_PI / 2
                 End If
      If (x_re < 0) And (x_im >= 0) Then
                 Arg = Atn(x_im / x_re) + M_PI
                 End If
      If (x_re < 0) And (x_im < 0) Then
                 Arg = Atn(x_im / x_re) - M_PI
                 End If
      If (x_re = 0) And (x_im = 0) Then
               '  MsgBox (" Xre=0 and Xim=0, Invalid result")
                 Arg = 0
                 End If

     End Function
     '36.  - 4. + 9.6568542i  - 4. + 4.i  - 4. + 1.6568542i  - 4.  - 4. - 1.6568542i
     
     '         column 7 to 8
     
      ' - 4. - 4.i  - 4. - 9.6568542i
     Sub GetIQfromR_fi(ByVal Num As Long, ByRef ModXk() As Double, ByRef argXk() As Double, _
     ByRef Xkre() As Double, ByRef Xkim() As Double, ByVal mode As Long)
       Dim i As Long
       Dim rad As Double
      
        rad = 1
         If mode = 1 Then
           rad = M_PI / 180
         End If
        
         For i = 1 To Num
          Xkre(i) = ModXk(i) * Cos(rad * argXk(i))
          Xkim(i) = ModXk(i) * Sin(rad * argXk(i))
         Next i
        
       End Sub

     


     Sub DFT(ByVal numdft As Long, ByRef dft_xnre() As Double, ByRef dft_xnim() As Double, ByRef dft_Xkre() As Double, ByRef dft_Xkim() As Double)

      Dim i, k, N As Long
     
     Dim angle As Double

     For k = 0 To numdft - 1
      dft_Xkre(k + 1) = 0
      dft_Xkim(k + 1) = 0
       For N = 0 To numdft - 1
        angle = -2 * M_PI * k * N / numdft
       dft_Xkre(k + 1) = dft_Xkre(k + 1) + dft_xnre(N + 1) * Cos(angle) - dft_xnim(N + 1) * Sin(angle)
       dft_Xkim(k + 1) = dft_Xkim(k + 1) + dft_xnre(N + 1) * Sin(angle) + dft_xnim(N + 1) * Cos(angle)
        Next N
      Next k

    End Sub
     
     Sub ScaleDFT(ByVal Num As Long, ByRef ScXkre() As Double, ByRef ScXkim() As Double)

    Dim i  As Long
      For i = 1 To Num
         ScXkre(i) = ScXkre(i) / Num
         ScXkim(i) = ScXkim(i) / Num
      Next i
     End Sub
     Sub ScaleIQidft(ByVal Num As Long, ByRef ScXk_re() As Double, ByRef ScXk_im() As Double)
       Dim i  As Long
       For i = 1 To Num
        ScXk_re(i) = ScXk_re(i) * Num
        ScXk_im(i) = ScXk_im(i) * Num
       Next i
      End Sub

     Sub IDFT(ByVal numidft As Long, ByRef idft_Xkre() As Double, ByRef idft_Xkim() As Double, ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)
     Dim k, N    As Long
      For N = 0 To numidft - 1
      Dim angle As Double
       'idft_item(n) = numidft - 1
       idft_xnre(N + 1) = 0
       idft_xnim(N + 1) = 0
        For k = 0 To numidft - 1
        angle = 2 * M_PI * k * N / numidft
        idft_xnre(N + 1) = idft_xnre(N + 1) + idft_Xkre(k + 1) * Cos(angle) - idft_Xkim(k + 1) * Sin(angle)
        idft_xnim(N + 1) = idft_xnim(N + 1) + idft_Xkre(k + 1) * Sin(angle) + idft_Xkim(k + 1) * Cos(angle)
        Next k
       
       idft_xnre(N + 1) = idft_xnre(N + 1) / numidft
       idft_xnim(N + 1) = idft_xnim(N + 1) / numidft

     Next N
     
      'idft_xnre(n), idft_xnim(n), numidft
     End Sub

     

     Sub Hilbert(ByRef xre() As Double, ByVal Numb As Long, ByRef Hilbertre() As Double, ByRef Hilbertim() As Double)
      'MATLAB algorithm
       Dim k, i, no2 As Long
       Dim h(1 To NMAX) As Double
     
      ' Dim dft_xnre(1 To NMAX)     As Double
       Dim dft_xnim(1 To NMAX)     As Double
       Dim dft_Xkre(1 To NMAX)    As Double
       Dim dft_Xkim(1 To NMAX)    As Double
       Dim idft_Xkre(1 To NMAX)   As Double
       Dim idft_Xkim(1 To NMAX)   As Double
       Dim idft_xnre(1 To NMAX)   As Double
       Dim idft_xnim(1 To NMAX)   As Double

      If (Numb = 0) Then
         Hilbertre(1) = 0
         Hilbertim(1) = 0
       MsgBox (" Hilbert() :   Parsing error ")
        End If


      no2 = Int(Numb / 2)
     
      
      
      
       For i = 1 To Numb
     
          dft_xnim(i) = 0
       Next i

       Call DFT(Numb, xre, dft_xnim, dft_Xkre, dft_Xkim)

       i = 0
         Do While i <= Numb

            h(i + 1) = 0
             i = i + 1
        Loop

      'if  N is even
     
        If (Numb > 0) And (2 * no2) = Numb Then
            h(1) = 1
            h(no2 + 1) = 1
            For k = 2 To no2
               h(k) = 2
            Next k

       'if N is odd

        ElseIf (Numb > 0) Then

             h(1) = 1
              For k = 2 To (0.5 * (Numb + 1))
                h(k) = 2
              Next k
          
         End If


         For i = 1 To Numb
          idft_Xkre(i) = dft_Xkre(i) * h(i)
          idft_Xkim(i) = dft_Xkim(i) * h(i)
         Next i
     
        
         Call IDFT(Numb, idft_Xkre, idft_Xkim, idft_xnre, idft_xnim)
       
       For i = 1 To Numb
       Hilbertre(i) = idft_xnre(i)
       Hilbertim(i) = idft_xnim(i)
       Next i
     
      End Sub

     Sub GetAnalyticSignal(ByRef xnre() As Double, ByVal Numb As Long, ByRef x_re() As Double, ByRef x_im() As Double)

      'I=x(t)
       'Q=jHilb(t)
       'Hilb=Hilbert()
       Dim i As Long
       Dim Hilbertre(1 To NMAX) As Double
       Dim Hilbertim(1 To NMAX) As Double
       Call Hilbert(xnre, Numb, Hilbertre, Hilbertim)
       For i = 1 To Numb
        x_re(i) = xnre(i)
        x_im(i) = Hilbertim(i)
       Next i
      End Sub

     Sub getenvelope(ByRef xnre() As Double, ByVal Numb As Long, ByRef envelope() As Double)
       Dim i As Long
       Dim Hilbertre(1 To NMAX) As Double
       Dim Hilbertim(1 To NMAX) As Double
       Call Hilbert(xnre, Numb, Hilbertre, Hilbertim)
       For i = 1 To Numb
       envelope(i) = ((Hilbertre(i) ^ 2) + (Hilbertim(i) ^ 2)) ^ 0.5
       Next i
      End Sub


     Sub GetSpectrumIQ(ByRef xn() As Double, ByVal Num As Long, ByRef Xkre() As Double, ByRef Xkim() As Double)
      Dim x_re(1 To NMAX) As Double
      Dim x_im(1 To NMAX) As Double

      Call GetAnalyticSignal(xn, Num, x_re, x_im)
      Call DFT(Num, x_re, x_im, Xkre, Xkim)
       Call ScaleDFT(Num, Xkre, Xkim)

     End Sub
     Sub GetR_fi(ByVal Numb As Long, ByRef Xkre() As Double, ByRef Xkim() As Double, _
     ByRef ModXk() As Double, ByRef argXk() As Double, ByVal mode As Long)
       Dim i As Long
       For i = 1 To Numb
        ModXk(i) = (Xkre(i) ^ 2 + Xkim(i) ^ 2) ^ 0.5
        argXk(i) = Arg(Xkre(i), Xkim(i))

      If mode = 1 Then
        argXk(i) = argXk(i) * 180 / M_PI
       End If

      Next i

     End Sub


    Sub GetSpectrum(ByRef xn() As Double, ByVal Num As Long, ByRef ModXk() As Double, ByRef argXk() As Double, ByVal mode As Long)
      Dim x_re(1 To NMAX) As Double
      Dim x_im(1 To NMAX) As Double
      Dim Xkre(1 To NMAX) As Double
      Dim Xkim(1 To NMAX) As Double
          Call GetAnalyticSignal(xn, Num, x_re, x_im)
          Call DFT(Num, x_re, x_im, Xkre, Xkim)
          Call ScaleDFT(Num, Xkre, Xkim)
          Call GetR_fi(Num, Xkre, Xkim, ModXk, argXk, mode)
     
      End Sub
     
     
      Sub GetSpectrumfromIQ(ByRef x_re() As Double, ByRef x_im() As Double, _
      ByVal Num As Long, ByRef ModXk() As Double, ByRef argXk() As Double, ByVal mode As Long)
      
      
      Dim Xkre(1 To NMAX) As Double
      Dim Xkim(1 To NMAX) As Double
          
          Call DFT(Num, x_re, x_im, Xkre, Xkim)
          Call ScaleDFT(Num, Xkre, Xkim)
          Call GetR_fi(Num, Xkre, Xkim, ModXk, argXk, mode)
     
      End Sub


       Sub SynthSignalfromIQ(ByVal Num As Long, ByRef Xkre() As Double, ByRef Xkim() As Double, _
       ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)
       
               
              'Call ScaleIQidft(Num, Xkre, Xkim)
              Call IDFT(Num, Xkre, Xkim, idft_xnre, idft_xnim)
              Call ScaleIQidft(Num, idft_xnre, idft_xnim)
       End Sub

       Sub SynthSignalfromRFi(ByRef ModXk() As Double, ByRef argXk() As Double, _
        ByRef xnre() As Double, ByRef xnim() As Double, ByVal Num As Long, ByVal mode As Long)
        Dim Xk_re(1 To NMAX) As Double
        Dim Xk_im(1 To NMAX) As Double
         Call GetIQfromR_fi(Num, ModXk, argXk, Xk_re, Xk_im, mode)
         Call SynthSignalfromIQ(Num, Xk_re, Xk_im, xnre, xnim)
       End Sub

    Sub printcomplex(ByRef xnre() As Double, ByRef xnim() As Double, _
     ByVal Num As Long, ByVal Col1 As Long, ByVal Col2 As Long)
         Dim i As Long
         For i = 1 To Num
              Cells(i + 1, Col1).Value = xnre(i)
              Cells(i + 1, Col2).Value = xnim(i)
            Next i
     End Sub

    Sub printarray(ByRef xnre() As Double, ByVal Num As Long, ByVal Col1 As Long)
         For i = 1 To Num
            Cells(i + 1, Col1).Value = xnre(i)
            
            Next i
     End Sub

     Function deltat(ByRef argF() As Double, ByVal f0 As Double, ByVal fs As Double, ByVal N As Long) As Double
      Dim fv(1 To NMAX + 1) As Double
      Dim N0 As Long
      Dim i As Long
      Call Getfassoc(fs, fv, N)
     N0 = 0
     For i = 1 To N
      If (fv(i) >= f0) And (fv(i + 1) <= f0) Then
       N0 = i
       End If
       Next i
     
     deltat = -(argF(N0) - argF(0)) / (2 * M_PI * fv(N0) - 2 * M_PI * fv(0))
     End Function


     Sub writetofile(ByRef xnre() As Double, ByRef xnim() As Double, _
     ByRef FilePath As String, ByVal N As Long)
      Dim FilePath As String
      Dim i As Long

     Open FilePath For Output As #1
     
         'Write #1,message
         'Write #1,"\n ; ; ; ;   "
         ' If  mode=0  Then
             'Write #1,"\n n ; Re(x) ; +j* Im(x)      ;   "
             '  End If
         
          'If  mode=1 Then
             ' Write #1,"\n n ; |X|   ; exp(j* Arg(x) );   "
             '  End If

        For i = 1 To N
             Write #1, i
             Write #1, xnre(i)
             Write #1, xnim(i)
         Next i
            
         Write #1, "\n ; ; ; ;  "
         Close #1
         MsgBox (" data saved to file")
       
     End Sub

     Sub TestHilbert()
    Dim fs As Double
     Dim fv(1 To NMAX) As Double
     Dim tv(1 To NMAX) As Double
     Dim Xi(1 To NMAX) As Double
     Dim Xq(1 To NMAX) As Double

     Dim Yi(1 To NMAX) As Double
     Dim Yq(1 To NMAX) As Double
     Dim Hi(1 To NMAX) As Double
     Dim Hq(1 To NMAX) As Double
     Dim Ii(1 To NMAX) As Double
     Dim Iq(1 To NMAX) As Double
     Dim modX(1 To NMAX) As Double
     Dim argX(1 To NMAX) As Double
     Dim env(1 To NMAX) As Double
     Dim modX1(1 To NMAX) As Double
     Dim argX1(1 To NMAX) As Double
     Dim Xi1(1 To NMAX) As Double
     Dim i As Long
     Dim N As Long
    Dim x_re(1 To NMAX) As Double
    Dim Yre(1 To NMAX) As Double
    Dim Yim(1 To NMAX) As Double
    N = 8
     Xi(1) = 1
     Xi(2) = 2
     Xi(3) = 3
     Xi(4) = 4
     Xi(5) = 5
     Xi(6) = 6
     Xi(7) = 7
     Xi(8) = 8


     Xq(1) = 0
     Xq(2) = 0
     Xq(3) = 0
     Xq(4) = 0
     Xq(5) = 0
     Xq(6) = 0
     Xq(7) = 0
     Xq(8) = 0
     fs = 1000
      Call printcomplex(Xi, Xq, N, 1, 2)
      Call Getfassoc(fs, fv, N)
      Call Gettassoc(fs, tv, N)
      Call printcomplex(tv, fv, N, 3, 4)
     
     Call DFT(N, Xi, Xq, Yi, Yq)
     
     Call printcomplex(Yi, Yq, N, 5, 6)
     
       Call IDFT(N, Yi, Yq, Ii, Iq)
        Call printcomplex(Ii, Iq, N, 7, 8)
       Call Hilbert(Xi, N, Hi, Hq)
        Call printcomplex(Hi, Hq, N, 9, 10)
        Call GetSpectrum(Xi, N, modX, argX, 1)
         Call printcomplex(modX, argX, N, 11, 12)
      
     
       Call getenvelope(Xi, N, env)
       Call printarray(env, N, 13)
     
     
      
       N = 16
     
     
       For i = 0 To N - 1
        argX(i + 1) = 0
        modX(i + 1) = 0
        Next i
      
       modX(2) = 1
       argX(2) = -45
       modX(6) = 0
       argX(6) = 0
      
      '  cos (90° – x) = sin x
      '  -sin ( x-90) = cos x
      '
     
      ' Hilbert(sin(x)) = -cos(t)
      ' Hilbert(Cos(x)) = sin(x)
     
     
      ' For i = 0 To n - 1
      '   Xi1(i + 1) = 1 * Sin(2 * M_PI * 15.625 * i / fs + 0) + 0.5 * Cos(2 * M_PI * 93.75 * i / fs + 0)
      
      '  Next i
      
       Call SynthSignalfromRFi(modX, argX, Xi, Xq, N, 1)
       Call printcomplex(Xi, Xq, N, 14, 15)
      
      
          Call Getfassoc(fs, fv, N)
          Call printarray(fv, N, 16)
         
         Call GetSpectrumfromIQ(Xi, Xq, N, modX1, argX1, 1)
         Call printcomplex(modX1, argX1, N, 17, 18)
         
         For i = 0 To N - 1
          Xi1(i + 1) = 0.5 * Cos(2 * M_PI * 250 * i / fs + 0) + 1 * Sin(2 * M_PI * 125 * i / fs + 0)
      
         Next i
          'Xi1(1) = 0
         
         Call GetSpectrum(Xi1, N, modX1, argX1, 1)
         Call printcomplex(modX1, argX1, N, 19, 20)
     
         Call GetSpectrumIQ(Xi1, N, Yre, Yim)
         Call printcomplex(Yre, Yim, N, 21, 22)
     
     End Sub

    Function GetNumF1(ByVal f1 As Double, ByVal Numb As Long, ByVal fs As Double) As Long

     GetNumF1 = 1 + Int(f1 * Numb / fs)

    End Function

    Function GetTHD_F(ByRef xn() As Double, ByVal Numb As Long, ByVal fs As Double, ByVal f1 As Double) As Double


    Dim numbf1 As Long
    Dim ModXk(1 To NMAX) As Double
    Dim argXk(1 To NMAX) As Double
    Dim Urms1 As Double
    Dim Uhharm   As Double
    Dim i As Long


    numbf1 = Int(f1 * Numb / fs) ' 2  , get associated vector of frequency, ORIGIN=0
    Call GetSpectrum(xn, Numb, ModXk, argXk, 0)
    Urms1 = ModXk(numbf1 + 1) / ((2) ^ (1 / 2)) ' / 1.41

    i = numbf1 * 2
    Uhharm = 0
    Do While (i <= Numb + 1)
     If (i + 1 <= Numb) Then
     Uhharm = Uhharm + ((ModXk(i + 1) ^ 2) / 2)    ' 2*2 +1, 3*2+1
     End If
    i = i + numbf1
    Loop

    Uhharm = (Uhharm) ^ (1 / 2)
    GetTHD_F = Uhharm * 100 / (Urms1)

    End Function


    Function GetTHD_R(ByRef xn() As Double, ByVal Numb As Long, ByVal fs As Double, ByVal f1 As Double) As Double

    Dim numbf1 As Long
    Dim ModXk(1 To NMAX) As Double
    Dim argXk(1 To NMAX) As Double
    Dim Urmssigah As Double
    Dim Uhharm   As Double
    Dim i As Long


     
    numbf1 = Int(f1 * Numb / fs)    ' 2 , get associated vector of frequency, ORIGIN=0
     
    Call GetSpectrum(xn, Numb, ModXk, argXk, 0)

    Urmssigah = 0
    i = numbf1

    Do While (i <= Numb)

    If (i + 1 <= Numb) Then
    Urmssigah = Urmssigah + (ModXk(i + 1) ^ 2) / 2
    End If
    i = i + numbf1
    Loop

    Urmssigah = (Urmssigah) ^ (0.5)

    Uhharm = 0


    ' 2*2=4

    i = numbf1 * 2 ' number of  bin of 2-nd harmonic , for example, 4
    Do While (i <= Numb)
                            ' ModXk(4+1=5)
    If (i + 1 <= Numb) Then
    Uhharm = Uhharm + (ModXk(i + 1) ^ 2) / 2
    End If
    i = i + numbf1          ' 4+2=6,next harmonic= ModXk(6+1=7)
    Loop
    Uhharm = (Uhharm) ^ (0.5)

    GetTHD_R = Uhharm * 100 / Urmssigah


    End Function

     Function GetRMS(ByRef xn() As Double, ByVal Numb As Long)
    Dim i As Long
    Dim Urms As Double
    Urms = 0
    For i = 1 To Numb
     Urms = Urms + xn(i) ^ 2
    Next i
    Urms = Urms / Numb
    GetRMS = (Urms) ^ (0.5)


    End Function


    Function GetTHD_N(ByRef xn() As Double, ByVal Numb As Long, ByVal fs As Double, ByVal f1 As Double)

    Dim numbf1 As Long
    Dim ModXk(1 To NMAX) As Double
    Dim argXk(1 To NMAX) As Double
    Dim Urmssig As Double
    Dim Uhharmn   As Double
    Dim i As Long

     

    Urmssig = 0
    For i = 1 To Numb
     Urmssig = Urmssig + xn(i) ^ 2
    Next i
    Urmssig = Urmssig / Numb
    Urmssig = (Urmssig) ^ (0.5)
    'Cells(5, 1).Value = Urmssig

     
    numbf1 = 1 + Int(f1 * Numb / fs)
    'Cells(6, 1).Value = numbf1

    Call GetSpectrum(xn, Numb, ModXk, argXk, 0)

    Uhharmn = 0
    i = 1

    Do While (i <= Numb)
    Cells(i + 1, 7) = ModXk(i)
    If (i <> numbf1) Then


    Uhharmn = Uhharmn + 0.5 * (ModXk(i) ^ 2)

    End If

    i = i + 1

    Loop
    Uhharmn = (Uhharmn) ^ (0.5)
    'Cells(13, 1) = Uhharmn
    GetTHD_N = Uhharmn * 100 / Urmssig


    End Function


    Sub TestTHDN()
    Dim modX(1 To NMAX) As Double
    Dim argX(1 To NMAX) As Double
    Dim Xi(1 To NMAX)   As Double
    Dim Xq(1 To NMAX)   As Double
    Dim fvect(1 To NMAX) As Double
    Dim fs As Double
    Dim f1 As Double
    Dim N As Long
    Dim i As Long

    N = Cells(1, 1).Value
    fs = Cells(2, 1).Value
    f1 = Cells(3, 1).Value
    Call Getfassoc(fs, fvect, N)

    For i = 1 To N

     Cells(i + 1, 2).Value = fvect(i)

    modX(i) = Cells(i + 1, 3).Value
    argX(i) = Cells(i + 1, 4).Value

    Next i


      Call SynthSignalfromRFi(modX, argX, Xi, Xq, N, 1)
      Call printcomplex(Xi, Xq, N, 5, 6)
     
     
     Cells(15, 1).Value = GetTHD_N(Xi, N, fs, f1)
     Cells(16, 1).Value = GetTHD_R(Xi, N, fs, f1)
     Cells(17, 1).Value = GetTHD_F(Xi, N, fs, f1)


    End Sub

    Friday, November 27, 2015 8:47 PM
  • 16 F F fi s(t) sq(t) Ui
    10000 0 0 0 1,1 0 1,08E-15 0 1
    1250 625 0 0 0,707107 0,807107 2,15E-15 1 2 1
    1250 1 0 -0,1 1 1 2 3 1 2
    0,710634 1875 0 0 -0,70711 0,607107 1,95E-15 3 4 3
    3 2500 0,1 0 -0,9 2,58E-15 0,1 4 5 2 4
    3125 0 0 -0,70711 -0,60711 2,3E-15 5 6 5
    3750 0 0 -0,1 -1 2,87E-15 6 7 3 6
       Кнопка 1      
    4375 0 0 0,707107 -0,80711 2,63E-15 7 8 7
    5000 0 0 1,1 -7,8E-15 3,17E-15 8 9 4 8
    f 5625 0 0 0,707107 0,807107 3,16E-15 9 10 9
    6250 0 0 -0,1 1 3,25E-15 10 11 5 10
    0,070711 6875 0 0 -0,70711 0,607107 3,25E-15 11 12 11
    7500 0 0 -0,9 7,4E-15 3,43E-15 12 13 12
    9,950372 8125 0 0 -0,70711 -0,60711 3,57E-15 13 14 6 13
    9,950372 8750 0 0 -0,1 -1 4,07E-15 14 15 14
    10 9375 0 0 0,707107 -0,80711 5,93E-15 15 16 7 15
    Friday, November 27, 2015 8:47 PM
  •  Parser for CF=Xpeak/Xrms, THD+N,THD-R,THD-F, Xmid, Xrms ,Q , CFP=Pmax/Pmid

    Option Explicit
     Const M_PI As Double = 3.14159265358979
     
      Const NMAX As Long = 131072

     Sub Getfassoc(ByVal fs As Double, ByRef fvect() As Double, ByVal Numb As Long)
       Dim k As Long
        For k = 0 To Numb - 1
         fvect(k + 1) = fs * k / Numb
        Next k
      End Sub
     
     
      Sub Gettassoc(ByVal fs As Double, ByRef tvect() As Double, ByVal Numb As Long)
       Dim N As Long
        For N = 0 To Numb - 1
         tvect(N + 1) = N / fs
        Next N
      End Sub


      Function Arg(ByVal x_re As Double, ByVal x_im As Double) As Double
       'Dim x_re,x_im As Double
       'Dim arg As Double
       Dim message As String
      If (x_re > 0) Then
                 Arg = Atn(x_im / x_re)
                 End If
      If (x_re = 0) And (x_im > 0) Then
                 Arg = M_PI / 2
                 End If
      If (x_re = 0) And (x_im < 0) Then
                 Arg = -M_PI / 2
                 End If
      If (x_re < 0) And (x_im >= 0) Then
                 Arg = Atn(x_im / x_re) + M_PI
                 End If
      If (x_re < 0) And (x_im < 0) Then
                 Arg = Atn(x_im / x_re) - M_PI
                 End If
      If (x_re = 0) And (x_im = 0) Then
               '  MsgBox (" Xre=0 and Xim=0, Invalid result")
                 Arg = 0
                 End If

     End Function
     '36.  - 4. + 9.6568542i  - 4. + 4.i  - 4. + 1.6568542i  - 4.  - 4. - 1.6568542i
     
     '         column 7 to 8
     
      ' - 4. - 4.i  - 4. - 9.6568542i
     Sub GetIQfromR_fi(ByVal Num As Long, ByRef ModXk() As Double, ByRef argXk() As Double, _
     ByRef Xkre() As Double, ByRef Xkim() As Double, ByVal mode As Long)
       Dim i As Long
       Dim rad As Double
      
        rad = 1
         If mode = 1 Then
           rad = M_PI / 180
         End If
        
         For i = 1 To Num
          Xkre(i) = ModXk(i) * Cos(rad * argXk(i))
          Xkim(i) = ModXk(i) * Sin(rad * argXk(i))
         Next i
        
       End Sub

     


     Sub DFT(ByVal numdft As Long, ByRef dft_xnre() As Double, ByRef dft_xnim() As Double, ByRef dft_Xkre() As Double, ByRef dft_Xkim() As Double)

      Dim i, k, N As Long
     
     Dim angle As Double

     For k = 0 To numdft - 1
      dft_Xkre(k + 1) = 0
      dft_Xkim(k + 1) = 0
       For N = 0 To numdft - 1
        angle = -2 * M_PI * k * N / numdft
       dft_Xkre(k + 1) = dft_Xkre(k + 1) + dft_xnre(N + 1) * Cos(angle) - dft_xnim(N + 1) * Sin(angle)
       dft_Xkim(k + 1) = dft_Xkim(k + 1) + dft_xnre(N + 1) * Sin(angle) + dft_xnim(N + 1) * Cos(angle)
        Next N
      Next k

    End Sub
     
     Sub ScaleDFT(ByVal Num As Long, ByRef ScXkre() As Double, ByRef ScXkim() As Double)

    Dim i  As Long
      For i = 1 To Num
         ScXkre(i) = ScXkre(i) / Num
         ScXkim(i) = ScXkim(i) / Num
      Next i
     End Sub
     Sub ScaleIQidft(ByVal Num As Long, ByRef ScXk_re() As Double, ByRef ScXk_im() As Double)
       Dim i  As Long
       For i = 1 To Num
        ScXk_re(i) = ScXk_re(i) * Num
        ScXk_im(i) = ScXk_im(i) * Num
       Next i
      End Sub

     Sub IDFT(ByVal numidft As Long, ByRef idft_Xkre() As Double, ByRef idft_Xkim() As Double, ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)
     Dim k, N    As Long
      For N = 0 To numidft - 1
      Dim angle As Double
       'idft_item(n) = numidft - 1
       idft_xnre(N + 1) = 0
       idft_xnim(N + 1) = 0
        For k = 0 To numidft - 1
        angle = 2 * M_PI * k * N / numidft
        idft_xnre(N + 1) = idft_xnre(N + 1) + idft_Xkre(k + 1) * Cos(angle) - idft_Xkim(k + 1) * Sin(angle)
        idft_xnim(N + 1) = idft_xnim(N + 1) + idft_Xkre(k + 1) * Sin(angle) + idft_Xkim(k + 1) * Cos(angle)
        Next k
       
       idft_xnre(N + 1) = idft_xnre(N + 1) / numidft
       idft_xnim(N + 1) = idft_xnim(N + 1) / numidft

     Next N
     
      'idft_xnre(n), idft_xnim(n), numidft
     End Sub

     

     Sub Hilbert(ByRef xre() As Double, ByVal Numb As Long, ByRef Hilbertre() As Double, ByRef Hilbertim() As Double)
      'MATLAB algorithm
       Dim k, i, no2 As Long
       Dim h(1 To NMAX) As Double
     
      ' Dim dft_xnre(1 To NMAX)     As Double
       Dim dft_xnim(1 To NMAX)     As Double
       Dim dft_Xkre(1 To NMAX)    As Double
       Dim dft_Xkim(1 To NMAX)    As Double
       Dim idft_Xkre(1 To NMAX)   As Double
       Dim idft_Xkim(1 To NMAX)   As Double
       Dim idft_xnre(1 To NMAX)   As Double
       Dim idft_xnim(1 To NMAX)   As Double

      If (Numb = 0) Then
         Hilbertre(1) = 0
         Hilbertim(1) = 0
       MsgBox (" Hilbert() :   Parsing error ")
        End If


      no2 = Int(Numb / 2)
     
      
      
      
       For i = 1 To Numb
     
          dft_xnim(i) = 0
       Next i

       Call DFT(Numb, xre, dft_xnim, dft_Xkre, dft_Xkim)

       i = 0
         Do While i <= Numb

            h(i + 1) = 0
             i = i + 1
        Loop

      'if  N is even
     
        If (Numb > 0) And (2 * no2) = Numb Then
            h(1) = 1
            h(no2 + 1) = 1
            For k = 2 To no2
               h(k) = 2
            Next k

       'if N is odd

        ElseIf (Numb > 0) Then

             h(1) = 1
              For k = 2 To (0.5 * (Numb + 1))
                h(k) = 2
              Next k
          
         End If


         For i = 1 To Numb
          idft_Xkre(i) = dft_Xkre(i) * h(i)
          idft_Xkim(i) = dft_Xkim(i) * h(i)
         Next i
     
        
         Call IDFT(Numb, idft_Xkre, idft_Xkim, idft_xnre, idft_xnim)
       
       For i = 1 To Numb
       Hilbertre(i) = idft_xnre(i)
       Hilbertim(i) = idft_xnim(i)
       Next i
     
      End Sub

     Sub GetAnalyticSignal(ByRef xnre() As Double, ByVal Numb As Long, ByRef x_re() As Double, ByRef x_im() As Double)

      'I=x(t)
       'Q=jHilb(t)
       'Hilb=Hilbert()
       Dim i As Long
       Dim Hilbertre(1 To NMAX) As Double
       Dim Hilbertim(1 To NMAX) As Double
       Call Hilbert(xnre, Numb, Hilbertre, Hilbertim)
       For i = 1 To Numb
        x_re(i) = xnre(i)
        x_im(i) = Hilbertim(i)
       Next i
      End Sub

     Sub getenvelope(ByRef xnre() As Double, ByVal Numb As Long, ByRef envelope() As Double)
       Dim i As Long
       Dim Hilbertre(1 To NMAX) As Double
       Dim Hilbertim(1 To NMAX) As Double
       Call Hilbert(xnre, Numb, Hilbertre, Hilbertim)
       For i = 1 To Numb
       envelope(i) = ((Hilbertre(i) ^ 2) + (Hilbertim(i) ^ 2)) ^ 0.5
       Next i
      End Sub


     Sub GetSpectrumIQ(ByRef xn() As Double, ByVal Num As Long, ByRef Xkre() As Double, ByRef Xkim() As Double)
      Dim x_re(1 To NMAX) As Double
      Dim x_im(1 To NMAX) As Double

      Call GetAnalyticSignal(xn, Num, x_re, x_im)
      Call DFT(Num, x_re, x_im, Xkre, Xkim)
       Call ScaleDFT(Num, Xkre, Xkim)

     End Sub
     Sub GetR_fi(ByVal Numb As Long, ByRef Xkre() As Double, ByRef Xkim() As Double, _
     ByRef ModXk() As Double, ByRef argXk() As Double, ByVal mode As Long)
       Dim i As Long
       For i = 1 To Numb
        ModXk(i) = (Xkre(i) ^ 2 + Xkim(i) ^ 2) ^ 0.5
        argXk(i) = Arg(Xkre(i), Xkim(i))

      If mode = 1 Then
        argXk(i) = argXk(i) * 180 / M_PI
       End If

      Next i

     End Sub


    Sub GetSpectrum(ByRef xn() As Double, ByVal Num As Long, ByRef ModXk() As Double, ByRef argXk() As Double, ByVal mode As Long)
      Dim x_re(1 To NMAX) As Double
      Dim x_im(1 To NMAX) As Double
      Dim Xkre(1 To NMAX) As Double
      Dim Xkim(1 To NMAX) As Double
          Call GetAnalyticSignal(xn, Num, x_re, x_im)
          Call DFT(Num, x_re, x_im, Xkre, Xkim)
          Call ScaleDFT(Num, Xkre, Xkim)
          Call GetR_fi(Num, Xkre, Xkim, ModXk, argXk, mode)
     
      End Sub
     
     
      Sub GetSpectrumfromIQ(ByRef x_re() As Double, ByRef x_im() As Double, _
      ByVal Num As Long, ByRef ModXk() As Double, ByRef argXk() As Double, ByVal mode As Long)
      
      
      Dim Xkre(1 To NMAX) As Double
      Dim Xkim(1 To NMAX) As Double
          
          Call DFT(Num, x_re, x_im, Xkre, Xkim)
          Call ScaleDFT(Num, Xkre, Xkim)
          Call GetR_fi(Num, Xkre, Xkim, ModXk, argXk, mode)
     
      End Sub


       Sub SynthSignalfromIQ(ByVal Num As Long, ByRef Xkre() As Double, ByRef Xkim() As Double, _
       ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)
       
               
              'Call ScaleIQidft(Num, Xkre, Xkim)
              Call IDFT(Num, Xkre, Xkim, idft_xnre, idft_xnim)
              Call ScaleIQidft(Num, idft_xnre, idft_xnim)
       End Sub

       Sub SynthSignalfromRFi(ByRef ModXk() As Double, ByRef argXk() As Double, _
        ByRef xnre() As Double, ByRef xnim() As Double, ByVal Num As Long, ByVal mode As Long)
        Dim Xk_re(1 To NMAX) As Double
        Dim Xk_im(1 To NMAX) As Double
         Call GetIQfromR_fi(Num, ModXk, argXk, Xk_re, Xk_im, mode)
         Call SynthSignalfromIQ(Num, Xk_re, Xk_im, xnre, xnim)
       End Sub

    Sub printcomplex(ByRef xnre() As Double, ByRef xnim() As Double, _
     ByVal Num As Long, ByVal Col1 As Long, ByVal Col2 As Long)
         Dim i As Long
         For i = 1 To Num
              Cells(i + 1, Col1).Value = xnre(i)
              Cells(i + 1, Col2).Value = xnim(i)
            Next i
     End Sub

    Sub printarray(ByRef xnre() As Double, ByVal Num As Long, ByVal Col1 As Long)
         For i = 1 To Num
            Cells(i + 1, Col1).Value = xnre(i)
            
            Next i
     End Sub

     Function deltat(ByRef argF() As Double, ByVal f0 As Double, ByVal fs As Double, ByVal N As Long) As Double
      Dim fv(1 To NMAX + 1) As Double
      Dim N0 As Long
      Dim i As Long
      Call Getfassoc(fs, fv, N)
     N0 = 0
     For i = 1 To N
      If (fv(i) >= f0) And (fv(i + 1) <= f0) Then
       N0 = i
       End If
       Next i
     
     deltat = -(argF(N0) - argF(0)) / (2 * M_PI * fv(N0) - 2 * M_PI * fv(0))
     End Function


     Sub writetofile(ByRef xnre() As Double, ByRef xnim() As Double, _
     ByRef FilePath As String, ByVal N As Long)
      Dim FilePath As String
      Dim i As Long

     Open FilePath For Output As #1
     
         'Write #1,message
         'Write #1,"\n ; ; ; ;   "
         ' If  mode=0  Then
             'Write #1,"\n n ; Re(x) ; +j* Im(x)      ;   "
             '  End If
         
          'If  mode=1 Then
             ' Write #1,"\n n ; |X|   ; exp(j* Arg(x) );   "
             '  End If

        For i = 1 To N
             Write #1, i
             Write #1, xnre(i)
             Write #1, xnim(i)
         Next i
            
         Write #1, "\n ; ; ; ;  "
         Close #1
         MsgBox (" data saved to file")
       
     End Sub

     Sub TestHilbert()
    Dim fs As Double
     Dim fv(1 To NMAX) As Double
     Dim tv(1 To NMAX) As Double
     Dim Xi(1 To NMAX) As Double
     Dim Xq(1 To NMAX) As Double

     Dim Yi(1 To NMAX) As Double
     Dim Yq(1 To NMAX) As Double
     Dim Hi(1 To NMAX) As Double
     Dim Hq(1 To NMAX) As Double
     Dim Ii(1 To NMAX) As Double
     Dim Iq(1 To NMAX) As Double
     Dim modX(1 To NMAX) As Double
     Dim argX(1 To NMAX) As Double
     Dim env(1 To NMAX) As Double
     Dim modX1(1 To NMAX) As Double
     Dim argX1(1 To NMAX) As Double
     Dim Xi1(1 To NMAX) As Double
     Dim i As Long
     Dim N As Long
    Dim x_re(1 To NMAX) As Double
    Dim Yre(1 To NMAX) As Double
    Dim Yim(1 To NMAX) As Double
    N = 8
     Xi(1) = 1
     Xi(2) = 2
     Xi(3) = 3
     Xi(4) = 4
     Xi(5) = 5
     Xi(6) = 6
     Xi(7) = 7
     Xi(8) = 8


     Xq(1) = 0
     Xq(2) = 0
     Xq(3) = 0
     Xq(4) = 0
     Xq(5) = 0
     Xq(6) = 0
     Xq(7) = 0
     Xq(8) = 0
     fs = 1000
      Call printcomplex(Xi, Xq, N, 1, 2)
      Call Getfassoc(fs, fv, N)
      Call Gettassoc(fs, tv, N)
      Call printcomplex(tv, fv, N, 3, 4)
     
     Call DFT(N, Xi, Xq, Yi, Yq)
     
     Call printcomplex(Yi, Yq, N, 5, 6)
     
       Call IDFT(N, Yi, Yq, Ii, Iq)
        Call printcomplex(Ii, Iq, N, 7, 8)
       Call Hilbert(Xi, N, Hi, Hq)
        Call printcomplex(Hi, Hq, N, 9, 10)
        Call GetSpectrum(Xi, N, modX, argX, 1)
         Call printcomplex(modX, argX, N, 11, 12)
      
     
       Call getenvelope(Xi, N, env)
       Call printarray(env, N, 13)
     
     
      
       N = 16
     
     
       For i = 0 To N - 1
        argX(i + 1) = 0
        modX(i + 1) = 0
        Next i
      
       modX(2) = 1
       argX(2) = -45
       modX(6) = 0
       argX(6) = 0
      
      '  cos (90° – x) = sin x
      '  -sin ( x-90) = cos x
      '
     
      ' Hilbert(sin(x)) = -cos(t)
      ' Hilbert(Cos(x)) = sin(x)
     
     
      ' For i = 0 To n - 1
      '   Xi1(i + 1) = 1 * Sin(2 * M_PI * 15.625 * i / fs + 0) + 0.5 * Cos(2 * M_PI * 93.75 * i / fs + 0)
      
      '  Next i
      
       Call SynthSignalfromRFi(modX, argX, Xi, Xq, N, 1)
       Call printcomplex(Xi, Xq, N, 14, 15)
      
      
          Call Getfassoc(fs, fv, N)
          Call printarray(fv, N, 16)
         
         Call GetSpectrumfromIQ(Xi, Xq, N, modX1, argX1, 1)
         Call printcomplex(modX1, argX1, N, 17, 18)
         
         For i = 0 To N - 1
          Xi1(i + 1) = 0.5 * Cos(2 * M_PI * 250 * i / fs + 0) + 1 * Sin(2 * M_PI * 125 * i / fs + 0)
      
         Next i
          'Xi1(1) = 0
         
         Call GetSpectrum(Xi1, N, modX1, argX1, 1)
         Call printcomplex(modX1, argX1, N, 19, 20)
     
         Call GetSpectrumIQ(Xi1, N, Yre, Yim)
         Call printcomplex(Yre, Yim, N, 21, 22)
     
     End Sub

    Function GetNumF1(ByVal f1 As Double, ByVal Numb As Long, ByVal fs As Double) As Long

     GetNumF1 = 1 + Int(f1 * Numb / fs)

    End Function

    Function GetTHD_F(ByRef xn() As Double, ByVal Numb As Long, ByVal fs As Double, ByVal f1 As Double) As Double


    Dim numbf1 As Long
    Dim ModXk(1 To NMAX) As Double
    Dim argXk(1 To NMAX) As Double
    Dim Urms1 As Double
    Dim Uhharm   As Double
    Dim i As Long


    numbf1 = Int(f1 * Numb / fs) ' 2  , get associated vector of frequency, ORIGIN=0
    Call GetSpectrum(xn, Numb, ModXk, argXk, 0)
    Urms1 = ModXk(numbf1 + 1) / ((2) ^ (1 / 2)) ' / 1.41

    i = numbf1 * 2
    Uhharm = 0
    Do While (i <= Numb + 1)
     If (i + 1 <= Numb) Then
     Uhharm = Uhharm + ((ModXk(i + 1) ^ 2) / 2)    ' 2*2 +1, 3*2+1
     End If
    i = i + numbf1
    Loop

    Uhharm = (Uhharm) ^ (1 / 2)
    GetTHD_F = Uhharm * 100 / (Urms1)

    End Function


    Function GetTHD_R(ByRef xn() As Double, ByVal Numb As Long, ByVal fs As Double, ByVal f1 As Double) As Double

    Dim numbf1 As Long
    Dim ModXk(1 To NMAX) As Double
    Dim argXk(1 To NMAX) As Double
    Dim Urmssigah As Double
    Dim Uhharm   As Double
    Dim i As Long


     
    numbf1 = Int(f1 * Numb / fs)    ' 2 , get associated vector of frequency, ORIGIN=0
     
    Call GetSpectrum(xn, Numb, ModXk, argXk, 0)

    Urmssigah = 0
    i = numbf1

    Do While (i <= Numb)

    If (i + 1 <= Numb) Then
    Urmssigah = Urmssigah + (ModXk(i + 1) ^ 2) / 2
    End If
    i = i + numbf1
    Loop

    Urmssigah = (Urmssigah) ^ (0.5)

    Uhharm = 0


    ' 2*2=4

    i = numbf1 * 2 ' number of  bin of 2-nd harmonic , for example, 4
    Do While (i <= Numb)
                            ' ModXk(4+1=5)
    If (i + 1 <= Numb) Then
    Uhharm = Uhharm + (ModXk(i + 1) ^ 2) / 2
    End If
    i = i + numbf1          ' 4+2=6,next harmonic= ModXk(6+1=7)
    Loop
    Uhharm = (Uhharm) ^ (0.5)

    GetTHD_R = Uhharm * 100 / Urmssigah


    End Function

     Function GetRMS(ByRef xn() As Double, ByVal Numb As Long) As Double
    Dim i As Long
    Dim Urms As Double
    Urms = 0
    For i = 1 To Numb
     Urms = Urms + xn(i) ^ 2
    Next i
    Urms = Urms / Numb
    GetRMS = (Urms) ^ (0.5)


    End Function

    Function GetMidValue(ByRef xn() As Double, ByVal Numb As Long) As Double
    Dim i As Long
    Dim Xmid As Double
    Xmid = 0
    For i = 1 To Numb
     Xmid = Xmid + xn(i)
    Next i
    GetMidValue = Xmid / Numb
     
    End Function

    Function GetPeak(ByRef xn() As Double, ByVal Numb As Long, ByVal mode As Integer) As Double
    Dim GetPeak_plus As Double
    Dim GetPeak_minus As Double
    Dim GetPeak_max As Double
    Dim i As Long
    GetPeak_plus = 0
    GetPeak_minus = 0
    GetPeak_max = 0

    For i = 1 To Numb
     
     If (xn(i) >= GetPeak_plus) Then
     GetPeak_plus = xn(i)
     Else
     GetPeak_plus = GetPeak_plus
     End If
     
     
     If (xn(i) <= GetPeak_minus) Then
     GetPeak_minus = xn(i)
     Else
     GetPeak_minus = GetPeak_minus
     End If
     
      If (Abs(xn(i)) >= Abs(GetPeak_max)) Then
     GetPeak_max = xn(i)
     Else
     GetPeak_max = GetPeak_max
     End If
     
    Next i

     If (mode = 1) Then
     GetPeak = GetPeak_plus
     End If
     If (mode = 2) Then
     GetPeak = GetPeak_minus
     
     End If
     
     If (mode = 0) Then ' ïèêîâîå ïî ìîäóëþ
     GetPeak = GetPeak_max
     End If
     
      If (mode = 4) Then ' ïèêîâîå ïî ìîäóëþ
     GetPeak = Abs(GetPeak_max)
     End If
      If (mode = 3) Then ' ðàçìàõ
     GetPeak = GetPeak_plus - GetPeak_minus
     End If


    End Function

     Function GetCF(ByRef xn() As Double, ByVal Numb As Long) As Double
     GetCF = Abs(GetPeak(xn, Numb, 0)) / GetRMS(xn, Numb)
     End Function
     
      Function GetPAR(ByRef xn() As Double, ByVal Numb As Long) As Double
     GetPAR = Abs(GetPeak(xn, Numb, 0)) / GetMidValue(xn, Numb)
     End Function

     Function GetQ(ByVal T As Double, ByVal tau As Double) As Double
     GetQ = T / tau
     End Function

    Function GetTHD_N(ByRef xn() As Double, ByVal Numb As Long, ByVal fs As Double, ByVal f1 As Double)

    Dim numbf1 As Long
    Dim ModXk(1 To NMAX) As Double
    Dim argXk(1 To NMAX) As Double
    Dim Urmssig As Double
    Dim Uhharmn   As Double
    Dim i As Long

    Urmssig = 0
    For i = 1 To Numb
     Urmssig = Urmssig + xn(i) ^ 2
    Next i
    Urmssig = Urmssig / Numb
    Urmssig = (Urmssig) ^ (0.5)
    'Cells(5, 1).Value = Urmssig

     
    numbf1 = 1 + Int(f1 * Numb / fs)
    'Cells(6, 1).Value = numbf1

    Call GetSpectrum(xn, Numb, ModXk, argXk, 0)

    Uhharmn = 0
    i = 1

    Do While (i <= Numb)
    Cells(i + 1, 4) = ModXk(i)
    If (i <> numbf1) Then


    Uhharmn = Uhharmn + 0.5 * (ModXk(i) ^ 2)

    End If

    i = i + 1

    Loop
    Uhharmn = (Uhharmn) ^ (0.5)
    'Cells(13, 1) = Uhharmn
    GetTHD_N = Uhharmn * 100 / Urmssig


    End Function


    Sub TestTHDN()
    Dim modX(1 To NMAX) As Double
    Dim argX(1 To NMAX) As Double
    Dim Xi(1 To NMAX)   As Double
    Dim Xq(1 To NMAX)   As Double
    Dim tvect(1 To NMAX) As Double
    Dim fvect(1 To NMAX) As Double
    Dim fs As Double
    Dim f1 As Double
    Dim N As Long
    Dim i As Long

    N = Cells(1, 1).Value
    fs = Cells(2, 1).Value
    f1 = Cells(3, 1).Value
    Call Gettassoc(fs, tvect, N)
    Call Getfassoc(fs, fvect, N)
    For i = 1 To N

     
    Cells(i + 1, 2).Value = tvect(i)
    Cells(i + 1, 5).Value = fvect(i)
    Xi(i) = Cells(i + 1, 3).Value


    Next i


     
     
     
     Cells(2, 8).Value = GetTHD_N(Xi, N, fs, f1)
     Cells(3, 8).Value = GetTHD_R(Xi, N, fs, f1)
     Cells(4, 8).Value = GetTHD_F(Xi, N, fs, f1)
     Cells(5, 8).Value = GetRMS(Xi, N)
     Cells(6, 8).Value = GetMidValue(Xi, N)
                       
      Cells(7, 8).Value = GetPeak(Xi, N, 0)
                         
      Cells(8, 8).Value = GetPeak(Xi, N, 1)
                         
      Cells(9, 8).Value = GetPeak(Xi, N, 2)
                         
      Cells(10, 8).Value = GetPeak(Xi, N, 3)
       Cells(11, 8).Value = GetCF(Xi, N)
       Cells(12, 8).Value = GetPAR(Xi, N)
    End Sub


    • Edited by USERPC01 Saturday, November 28, 2015 8:56 PM
    Friday, November 27, 2015 10:58 PM
  • 32 t s(t) Ui f
    192000 0 0 0,03125 0 THD+N 95,19716    Кнопка 1      
    6000 5,21E-06 0 0,0625 6000 THD-R 96,66572
    1,04E-05 0 0,0625 12000 THD-F 377,4917
    1,56E-05 0 0,0625 18000 Vrms 0,176777
    2,08E-05 0 0,0625 24000 Vmid 0,03125
    2,6E-05 0 0,0625 30000 |Vmax| 1
    3,13E-05 0 0,0625 36000 Vpeak+ 1
    3,65E-05 0 0,0625 42000 Vpeak- 0
    4,17E-05 0 0,0625 48000 размах 1
    4,69E-05 0 0,0625 54000 CF 5,656854
    5,21E-05 0 0,0625 60000 CFP 32
    5,73E-05 0 0,0625 66000
    6,25E-05 1 0,0625 72000
    6,77E-05 0 0,0625 78000
    7,29E-05 0 0,0625 84000
    7,81E-05 0 0,0625 90000
    8,33E-05 0 0,03125 96000
    8,85E-05 0 1,92E-16 102000
    9,38E-05 0 2,39E-16 108000
    9,9E-05 0 1,39E-16 114000
    0,000104 0 1,08E-16 120000
    0,000109 0 1,41E-16 126000
    0,000115 0 6,17E-17 132000
    0,00012 0 8,12E-17 138000
    0,000125 0 2,05E-16 144000
    0,00013 0 8,29E-17 150000
    0,000135 0 1,6E-16 156000
    0,000141 0 9,21E-17 162000
    0,000146 0 2,26E-16 168000
    0,000151 0 1,49E-16 174000
    0,000156 0 1,64E-16 180000
    0,000161 0 8,76E-16 186000
    Friday, November 27, 2015 11:00 PM
  • Hello,

    The 'Academic Project Program' forum is for posts related to Microsoft's Academic Project Program.

    Are you trying to ask questions about your code? if so, please ask in the Visual C++ Language forum on MSDN.

    If you just want to share your code, please use https://code.msdn.microsoft.com/

    Karl


    When you see answers and helpful posts, please click Vote As Helpful, Propose As Answer, and/or Mark As Answer.
    My Blog: Unlock PowerShell
    My Book: Windows PowerShell 2.0 Bible
    My E-mail: -join('6D73646E5F6B61726C406F75746C6F6F6B2E636F6D'-split'(?<=\G.{2})'|%{if($_){[char][int]"0x$_"}})

    Monday, November 30, 2015 5:23 PM
    Moderator