VB 5/6-Tipp 0716: Faltung und Filterung von Funktionen/Signalen
von Arne Elster
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: | Verwendete API-Aufrufe: keine | Download: |
'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-Version | Win32s | Win95 | Win98 | WinME | WinNT4 | Win2000 | WinXP |
VB4 | |||||||
VB5 | |||||||
VB6 |
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.