Die Community zu .NET und Classic VB.
Menü

VB 5/6-Tipp 0716: Faltung und Filterung von Funktionen/Signalen

 von 

Beschreibung 

Mit diesem Modul ist es möglich, höhere oder tiefere Frequenzbänder einer Funktion oder eines Signals zu filtern, also abzuschwächen.
Dabei kommt ein Windowed-Sinc-FIR-Filter zum Einsatz, dessen Qualität und Cutoff bestimmt werden können.
Das Filter wird mit der Funktion/dem Signal gefaltet, wobei hier eine Overlap-Technik zum Einsatz kommt, um fortlaufende Faltung zu ermöglichen.

Schwierigkeitsgrad:

Schwierigkeitsgrad 2

Verwendete API-Aufrufe:

keine

Download:

Download des Beispielprojektes [6,62 KB]

'Dieser Quellcode stammt von http://www.activevb.de
'und kann frei verwendet werden. Für eventuelle Schäden
'wird nicht gehaftet.

'Um Fehler oder Fragen zu klären, nutzen Sie bitte unser Forum.
'Ansonsten viel Spaß und Erfolg mit diesem Source!

'------------- Anfang Projektdatei Projekt1.vbp -------------
'--------- Anfang Formular "Form1" alias Form1.frm  ---------
' Steuerelement: Rahmensteuerelement "Frame3"
' Steuerelement: Bildfeld-Steuerelement "picHighPass" auf Frame3
' Steuerelement: Bildfeld-Steuerelement "picHPKernel" auf Frame3
' Steuerelement: Bildfeld-Steuerelement "picSpecHP" auf Frame3
' Steuerelement: Bildfeld-Steuerelement "picSpecHPK" auf Frame3
' Steuerelement: Beschriftungsfeld "lblKernel2" auf Frame3
' Steuerelement: Beschriftungsfeld "lblFreqSpec" (Index von 0 bis 2) auf Frame3
' Steuerelement: Rahmensteuerelement "Frame2"
' Steuerelement: Bildfeld-Steuerelement "picLowPass" auf Frame2
' Steuerelement: Bildfeld-Steuerelement "picLPKernel" auf Frame2
' Steuerelement: Bildfeld-Steuerelement "picSpecLP" auf Frame2
' Steuerelement: Bildfeld-Steuerelement "picSpecLPK" auf Frame2
' Steuerelement: Beschriftungsfeld "lblKernel" auf Frame2
' Steuerelement: Rahmensteuerelement "Frame1"
' Steuerelement: Schaltfläche "cmdRectangular" auf Frame1
' Steuerelement: Schaltfläche "cmdSinusoids" auf Frame1
' Steuerelement: Schaltfläche "cmdSawtooth" auf Frame1
' Steuerelement: Bildfeld-Steuerelement "picStd" auf Frame1
' Steuerelement: Kontrollkästchen-Steuerelement "chkNoise" auf Frame1
' Steuerelement: Bildfeld-Steuerelement "picSpecInp" auf Frame1
' Steuerelement: Vertikale Scrollbar "scrlTaps"
' Steuerelement: Vertikale Scrollbar "scrlFactor"
' Steuerelement: Beschriftungsfeld "lblTaps"
' Steuerelement: Beschriftungsfeld "lblFactor"
Option Explicit

Private Const SAMPLES           As Long = 512   ' Potenz von 2!

Private m_sngInp(SAMPLES - 1)   As Single

' Sägezahn Funktion, an der man mit einem Tiefpass
' die Entstehung der entsprechenden Fourier-Reihe
' betrachten kann
Private Function InputSawtooth(ByVal t As Single) As Single
    Const a As Single = 150
    
    InputSawtooth = 1.9 * (t / a - Fix(t / a + 0.5))
End Function

' Einfache Überlagerung von verschiedenen Sinusoiden
Private Function InputSinusoids(ByVal t As Single) As Single
    InputSinusoids = (Sin(2 * t) + Sin(0.3 * t) + Sin(0.1 * t)) * 0.5
End Function

Private Function InputRectangular(ByVal t As Single) As Single
    Const a As Single = 20
    Static x As Single
    
    If x = 0 Then x = 0.8
    
    If t Mod a = 0 Then
        x = -x
    End If
    
    InputRectangular = x
End Function

' Funktion/Signal skaliert in PictureBox plotten
Private Sub Plot(pb As PictureBox, sngSamples() As Single, _
    ByVal center As Boolean, ByVal normalize As Boolean)
    
    Dim dy      As Long, dy2        As Long
    Dim n       As Long, i          As Long
    Dim sngMax  As Single, sngVal   As Single
    Dim yL      As Single, yN       As Single
    Dim x       As Single, K        As Single
    Dim st      As Single
    
    dy = pb.ScaleHeight - 1
    dy2 = dy \ 2
    
    n = UBound(sngSamples) + 1
    st = n / pb.ScaleWidth
    
    If normalize Then
        For i = 0 To n - 1
            If Abs(sngSamples(i)) > sngMax Then sngMax = Abs(sngSamples(i))
        Next
    End If
    
    If sngMax = 0 Then sngMax = 1
    sngVal = sngSamples(0) / sngMax
    
    If center Then
        yL = -sngVal * dy2 + dy2
        
        pb.ForeColor = RGB(160, 160, 180)
        pb.Line (0, dy2)-(pb.ScaleWidth, dy2)
        pb.ForeColor = vbBlack
    Else
        yL = dy - sngVal * dy
    End If
    
    K = K + st
    
    Do
        sngVal = sngSamples(Fix(K)) / sngMax
        
        If center Then
            yN = -sngVal * dy2 + dy2
        Else
            yN = dy - sngVal * dy
        End If
        
        pb.Line (x, yL)-(x + 1, yN)
        yL = yN
        x = x + 1
        K = K + st
    Loop While K < n
End Sub

' Signal filtern und Ergebnis plotten
Private Sub Filter(ByVal lngTaps As Long, ByVal sngFactor As Single)
    Dim sngLP()         As Single
    Dim sngHP()         As Single
    Dim udtLP           As FilterKernel
    Dim udtHP           As FilterKernel
    
    udtLP = CreateFilter(FilterLowpass, lngTaps, sngFactor)
    udtHP = CreateFilter(FilterHighpass, lngTaps, sngFactor)
    
    sngLP = m_sngInp
    FilterProcess sngLP, udtLP
    picLowPass.Cls:     Plot picLowPass, sngLP, True, False
    picLPKernel.Cls:    Plot picLPKernel, udtLP.kernel, True, False
    picSpecLP.Cls:      DisplayFT picSpecLP, sngLP, False
    picSpecLPK.Cls:     DisplayFT picSpecLPK, udtLP.kernel, True
    
    sngHP = m_sngInp
    FilterProcess sngHP, udtHP
    picHighPass.Cls:    Plot picHighPass, sngHP, True, False
    picHPKernel.Cls:    Plot picHPKernel, udtHP.kernel, True, False
    picSpecHP.Cls:      DisplayFT picSpecHP, sngHP, False
    picSpecHPK.Cls:     DisplayFT picSpecHPK, udtHP.kernel, True
End Sub

' White Noise zum Signal dazurechnen
Private Sub chkNoise_Click()
    Dim i As Long
    
    If chkNoise.Value = 1 Then
        For i = 0 To SAMPLES - 1
            m_sngInp(i) = (m_sngInp(i) + Rnd() * 0.5) / 1.5
        Next
    Else
        For i = 0 To SAMPLES - 1
            m_sngInp(i) = InputSawtooth(i)
        Next
    End If

    UpdateInpDisplay
    UpdateDisplay
End Sub

Private Sub cmdRectangular_Click()
    Dim i As Long

    For i = 1 To SAMPLES - 1
        m_sngInp(i) = InputRectangular(i)
    Next
    
    UpdateInpDisplay
    UpdateDisplay
End Sub

Private Sub cmdSawtooth_Click()
    Dim i As Long
    
    For i = 0 To SAMPLES - 1
        m_sngInp(i) = InputSawtooth(i)
    Next
    
    UpdateInpDisplay
    UpdateDisplay
End Sub

Private Sub cmdSinusoids_Click()
    Dim i As Long
    
    For i = 0 To SAMPLES - 1
        m_sngInp(i) = InputSinusoids(i)
    Next
    
    UpdateInpDisplay
    UpdateDisplay
End Sub

Private Sub Form_Load()
    Dim i   As Long
    
    For i = 0 To SAMPLES - 1
        m_sngInp(i) = InputSawtooth(i)
    Next
    
    UpdateInpDisplay
    UpdateDisplay
End Sub

' Fourier Transformation von Daten, anschließend plotten
Private Sub DisplayFT(pb As PictureBox, data() As Single, ByVal normalize As Boolean)
    Dim sngReaInp(SAMPLES - 1)      As Single
    Dim sngRealOut(SAMPLES - 1)     As Single
    Dim sngImagOut(SAMPLES - 1)     As Single
    Dim sngCompOut(SAMPLES / 2 - 1) As Single
    Dim i                           As Long
    
    For i = 0 To UBound(data)
        sngReaInp(i) = data(i)
    Next
    
    RealFFT SAMPLES, sngReaInp, sngRealOut, sngImagOut
    
    For i = 0 To SAMPLES / 2 - 1
        sngCompOut(i) = Sqr(sngRealOut(i) * sngRealOut(i) + _
            sngImagOut(i) * sngImagOut(i)) / (SAMPLES / 2)
    Next
    
    Plot pb, sngCompOut, False, normalize
End Sub

Private Sub UpdateInpDisplay()
    picStd.Cls
    picSpecInp.Cls
    Plot picStd, m_sngInp, True, False
    DisplayFT picSpecInp, m_sngInp, False
End Sub

Private Sub UpdateDisplay()
    Dim sngFactor   As Single
    Dim lngTaps     As Long
    
    sngFactor = scrlFactor.Value / scrlFactor.Max * 0.5
    lngTaps = scrlTaps.Value
    
    lblFactor.Caption = "Faktor (" & Format(sngFactor, "0.00") & "):"
    lblTaps.Caption = "Taps (" & lngTaps & "):"
    Filter lngTaps, sngFactor
End Sub

Private Sub scrlFactor_Change()
    UpdateDisplay
End Sub

Private Sub scrlFactor_Scroll()
    UpdateDisplay
End Sub

Private Sub scrlTaps_Change()
    UpdateDisplay
End Sub

Private Sub scrlTaps_Scroll()
    UpdateDisplay
End Sub
'---------- Ende Formular "Form1" alias Form1.frm  ----------
'---- Anfang Modul "SignalFilter" alias SignalFilter.bas ----
Option Explicit


' Windowed Sinc Filter und Convolution
' nach http://www.dspguide.com/
'
' factor:   Wert von 0.0 bis 0.5, der den Cutoff angibt.
'           Bsp.: Audiosignal gesampelt mit 22050 Hz,
'                 Cutoff Frequenz soll bei 220.5 Hz liegen,
'                 dann wäre factor = 220.5/22050 = 0.01
'
' Cutoff:   Schwelle, ab der das Filter beginnt zu wirken
'
' Taps:     Qualität des Filters, je mehr Taps, desto steiler
'           das Rolloff Band nach dem Cutoff, aber je mehr Taps,
'           desto größer das Delay, und desto mehr Rechenarbeit


Public Type FilterKernel
    kernel()        As Single
    olap()          As Single
    taps            As Long
End Type

Public Enum FilterType
    FilterHighpass = 0
    FilterLowpass
End Enum

Private Const PI    As Single = 3.14159265358979
Private Const PI2   As Single = PI * 2

' Windows Sinc FIR Filter Kern berechnen
' nach http://www.dspguide.com/
Public Function CreateFilter(ByVal ftp As FilterType, ByVal taps As Long, _
    ByVal factor As Single) As FilterKernel
    
    Dim omega       As Single
    Dim n           As Single
    Dim sum         As Single
    Dim m           As Long
    Dim i           As Long

    omega = PI2 * factor
    m = taps / 2

    With CreateFilter
        ReDim .kernel(taps - 1) As Single
        ReDim .olap(taps - 1) As Single
        .taps = taps
    
        ' Sinc Funktion
        For i = 0 To taps - 1
            If i - m = 0 Then
                .kernel(i) = omega
            Else
                n = i - m
                .kernel(i) = Sin(omega * n) / n
            End If
            
            ' Glockenfenster um Ripple zu minimieren
            .kernel(i) = .kernel(i) * HammingWindow(i, taps)
            sum = sum + .kernel(i)
        Next
        
        If sum = 0 Then sum = 1
        
        For i = 0 To taps - 1
            ' Kern normalisieren
            .kernel(i) = .kernel(i) / sum
        Next
        
        If ftp = FilterHighpass Then
            ' Spektrale Inversion, um Lowpass in Highpass umzuwandeln
            For i = 0 To taps - 1
                .kernel(i) = -.kernel(i)
            Next
            .kernel(m) = .kernel(m) + 1
        End If
        
        '' durch Convolution des Kerns mit sich selbst
        '' kann das Stopband nochmal verstärkt werden,
        '' allerdings ergeben sich daraus auch mehr Taps
        ' Dim k2() As Single
        ' k2 = .kernel
        ' .kernel = Convolve(k2, .kernel)
        ' ReDim .olap(UBound(.kernel)) As Single
    End With
End Function

' Leert Overlap von Filter
Public Sub ResetFilter(kernel As FilterKernel)
    Dim i           As Long
    
    For i = 0 To kernel.taps - 1
        kernel.olap(i) = 0
    Next
End Sub

' Signalfaltung mit Filterkernel, Overlap-Add Methode
Public Sub FilterProcess(sngValues() As Single, kernel As FilterKernel)
    Dim sngOut()    As Single
    Dim i           As Long
    Dim n           As Long
    
    n = UBound(sngValues) + 1
    
    If n >= kernel.taps - 1 Then
        sngOut = Convolve(sngValues, kernel.kernel)

        For i = 0 To kernel.taps - 1
            sngOut(i) = sngOut(i) + kernel.olap(i)
            kernel.olap(i) = sngOut(n + i)
        Next

        For i = 0 To n - 1
            sngValues(i) = sngOut(i)
        Next
    End If
End Sub

' Faltung nach Input Side Methode
Private Function Convolve(a() As Single, B() As Single) As Single()
    Dim c() As Single
    Dim i   As Long
    Dim j   As Long
    
    ReDim c(UBound(a) + UBound(B) + 1) As Single
    
    For i = 0 To UBound(a)
        For j = 0 To UBound(B)
            c(i + j) = c(i + j) + a(i) * B(j)
        Next
    Next
    
    Convolve = c
End Function

' http://www.dspguide.com/
Private Function HammingWindow(ByVal i As Single, ByVal n As Single) As Single
    HammingWindow = 0.54 - 0.46 * Cos(PI2 * i / n)
End Function
'----- Ende Modul "SignalFilter" alias SignalFilter.bas -----
'------------- Anfang Modul "FFT" alias FFT.bas -------------
Option Explicit

' Radix-2 FFT von Murphy McCauley

Private Const PI    As Double = 3.14159265358979
Private Const PI2   As Double = PI * 2

Private m_lngP2(16) As Long

Private Function ReverseBits(ByVal Index As Long, NumBits As Long) As Long
    Dim i As Long, Rev As Long

    For i = 0& To NumBits - 1&
        Rev = (Rev * 2&) Or (Index And 1&)
        Index = Index \ 2&
    Next

    ReverseBits = Rev
End Function

Public Sub RealFFT( _
    ByVal NumSamples As Long, _
    RealIn() As Single, _
    RealOut() As Single, ImagOut() As Single, _
    Optional InverseTransform As Boolean = False _
)

    Dim AngleNumerator                  As Double

    Dim NumBits                         As Long

    Dim i           As Long, j          As Long
    Dim K           As Long, n          As Long

    Dim BlockSize   As Long, BlockEnd   As Long

    Dim DeltaAngle  As Single, DeltaAr  As Single
    Dim Alpha       As Single, Beta     As Single

    Dim TR          As Single, TI       As Single
    Dim AR          As Single, AI       As Single

    If m_lngP2(0) = 0 Then
        For i = 0 To 16
            m_lngP2(i) = 2 ^ i
        Next
    End If

    If InverseTransform Then
        AngleNumerator = -PI2
    Else
        AngleNumerator = PI2
    End If

    For i = 0 To 16
        If (NumSamples And m_lngP2(i)) <> 0 Then
            NumBits = i
            Exit For
        End If
    Next

    For i = 0 To (NumSamples - 1)
        j = ReverseBits(i, NumBits)
        RealOut(j) = RealIn(i)
        ImagOut(j) = 0
    Next

    BlockEnd = 1
    BlockSize = 2

    Do While BlockSize <= NumSamples
        DeltaAngle = AngleNumerator / BlockSize
        Alpha = Sin(0.5 * DeltaAngle)
        Alpha = 2# * Alpha * Alpha
        Beta = Sin(DeltaAngle)

        i = 0
        Do While i < NumSamples
            AR = 1#
            AI = 0#
            
            j = i
            For n = 0 To BlockEnd - 1
                K = j + BlockEnd
                TR = AR * RealOut(K) - AI * ImagOut(K)
                TI = AI * RealOut(K) + AR * ImagOut(K)
                RealOut(K) = RealOut(j) - TR
                ImagOut(K) = ImagOut(j) - TI
                RealOut(j) = RealOut(j) + TR
                ImagOut(j) = ImagOut(j) + TI
                DeltaAr = Alpha * AR + Beta * AI
                AI = AI - (Alpha * AI - Beta * AR)
                AR = AR - DeltaAr
                j = j + 1
            Next

            i = i + BlockSize

        Loop

        BlockEnd = BlockSize
        BlockSize = BlockSize * 2
    Loop

    If InverseTransform Then
        For i = 0 To NumSamples - 1
            RealOut(i) = RealOut(i) / NumSamples
            ImagOut(i) = ImagOut(i) / NumSamples
        Next
    End If
End Sub
'-------------- Ende Modul "FFT" alias FFT.bas --------------
'-------------- Ende Projektdatei Projekt1.vbp --------------

Tipp-Kompatibilität:

Windows/VB-VersionWin32sWin95Win98WinMEWinNT4Win2000WinXP
VB4
VB5
VB6

Hat dieser Tipp auf Ihrem Betriebsystem und mit Ihrer VB-Version funktioniert?

Ja, funktioniert!

Nein, funktioniert nicht bei mir!

VB-Version:

Windows-Version:

Ihre Meinung  

Falls Sie Fragen zu diesem Artikel haben oder Ihre Erfahrung mit anderen Nutzern austauschen möchten, dann teilen Sie uns diese bitte in einem der unten vorhandenen Themen oder über einen neuen Beitrag mit. Hierzu können sie einfach einen Beitrag in einem zum Thema passenden Forum anlegen, welcher automatisch mit dieser Seite verknüpft wird.