Ich möcht in diesem Beitrag zeigen, dass man Signalerarbeitung auch ohne viel Theorie und Mathematik hinbekommt und hab dafür eine Bibliothek geschrieben, die klein, supereinfach zu verwenden und vorallem auf einem ESP mit Micropython ihren Dienst tuen kann.

Signalverarbeitung war eines meiner Hauptfächer im Studium und ich kann mich noch an einen kompletten Leitzordner voller Theorie und Mathematik dazu erinnern. Sucht man im Internet nach z.B. einer FFT oder einen Filter in Python so findet man Programme die Bibliotheken verwenden, die es leider für Micropython nicht gibt und für die Ausführung auf einem Mikrocontroller nicht geeignet sind.

Eigentlich ist Signalverarbeitung recht einfach. Alles fing damit an, dass ein Französischer Mathematiker Jean Baptiste Joseph Fourier schon so um 1800 rum herausgefunden hat, dass sich jedes periodische Signal aus einer Überlagerung von Sinusschwingungen zusammensetzt.

Jede Sinusschwingung wird im Spektrum des Signals als Strich dargestellt. Hat man z.B. ein andersförmiges Signal, wie z.B. ein Rechteck, so besteht das Spektrum aus mehreren Strichen unterschiedlicher Frequenz und Amplitude.

Aber zurück zu unseren Sinusschwingungen. Ich hab ein kleines Programm geschrieben, das zwei Sinusschwingungen generiert und das zugehörige Spektrum anzeigt:

#####################################################################
# Testprogramm für b2c_sft : box2code - Signal Filter Toolbox
# Stefan am 2024-01-05
#####################################################################

######## Configure here:

sampleFreq = 2000
fftLen = 512 # must be power of 2
applyNotch = True
applyWinFilterBandStop = False 

########################


import matplotlib.pyplot as plt
import b2c_sft
import math

data1 = b2c_sft.generateSinus(fftLen, 100, 250, sampleFreq)
data2 = b2c_sft.generateSinus(fftLen, 50, 500, sampleFreq)

data = b2c_sft.mixSignals(data1,data2)
if applyNotch:
    notch = b2c_sft.Filter(b2c_sft.notchZeroes(500,sampleFreq), b2c_sft.notchPoles(500,sampleFreq))
    data = notch.filterIterative(data)

zeroes1 = [
    0.98811787286546693000, # comes from WinFilter
    0.00000000000000263287,
    2.96435361859640080000,
    0.00000000000000526575,
    2.96435361859640080000,
    0.00000000000000263287,
    0.98811787286546693000
]
poles1 = [
    1.00000000000000000000, # comes from WinFilter
    0.00000000013946281840,
    2.98743365009436030000,
    0.00000000021746471926,
    2.97494613251504750000,
    0.00000000001021009398,
    0.98751223597926108000
]

if applyWinFilterBandStop:
    bandStopFromWinFilter = b2c_sft.Filter(zeroes1, poles1)
    data = bandStopFromWinFilter.filterIterative(data)

plt.subplot(2, 1, 1) 
plt.plot(data)
plt.xlabel('Sample')
plt.ylabel('Amplitude')

b2c_sft.fft(data)

plt.subplot(2, 1, 2) 
plt.suptitle('Zeit- und Frequenzsignal')
plt.plot(b2c_sft.frequencyAxis(data, sampleFreq), data)
plt.xlabel('Frequenz in Hz')
plt.ylabel('Amplitude')

plt.show()

Ausschauen tut das ganze dann so (hab Thonny verwendet mit dem Lokalen Python-Interpreter. Die mathplotlib musst du dir vielleicht mit pip erst installieren – Der Code für die Bibliothek b2c_sft kommt etwas weiter unten):

Bei Bedarf kann man zwei Filter aktivieren. Der eine ist ein NotchFilter, der gezielt eine Sinusschwingung aus dem Spektrum herausfiltert:

Man sieht dass die Schwingung bei 500 Hz jetzt weg ist

Du hast hier bestimmt schon ein paar Fragen:

• Die Samplefrequenz ist doch 2000 – warum geht das Spektrum nur bis 1000

Zwei andere berühmte Mathematiker haben im Nyquist-Shannon-Abtasttheorem bestimmt, dass die Abtastfrequenz mindestens doppelt so hoch sein muss wie die höchste im Signal vorkommende Frequenz. D.h. der schnellste messbare Sinus darf höchstens halb so schnell sein wie die Samplefrequenz. Tatsächlich ist eine normale FFT um den Mittelpunkt spiegelsymmetrisch und enthält auch „negative Frequenzen“ aber weil ich diese Information in der Praxis für unnötig halte halbiere ich in der Lib die FFT nach der Berechnung. Dadurch gibt es keine Verwirrung.

• Das ist doch eine Filterbibliothek – Es gibt nur einen Notch-Filter – wie mach ich andere Filter wie Tief- oder Hochpass ?

Der Filter in der Bibliothek nimmt Listen von Polen und Nullstellen als Parameter (Eigentlich die Koeffizienten eines Polynoms mit entsprechenden Nullstellen, aber ich hab ja versprochen keine Mathematik 🙂 ). Das genügt um jeden erdenklichen Filter damit zu bauen. Aber wie kommt man an die Werte für Pole und Nullstellen ?

Dazu möchte ich dir WinFilter empfehlen:

Wenn du auf „C-Code Generieren“ klickst kommt als Ergebnis folgendes:

    float ACoef[NCoef+1] = {
        0.98811787286546693000,
        0.00000000000000263287,
        2.96435361859640080000,
        0.00000000000000526575,
        2.96435361859640080000,
        0.00000000000000263287,
        0.98811787286546693000
    };

    float BCoef[NCoef+1] = {
        1.00000000000000000000,
        0.00000000013946281840,
        2.98743365009436030000,
        0.00000000021746471926,
        2.97494613251504750000,
        0.00000000001021009398,
        0.98751223597926108000
    };

Hinweis: beim FIR-Filter gibt es nur einen Array – das sind dann deine zeroes – Übergeb dann für die pole einfach eine leere Liste.

Kopiere dir das raus und packe es, Python konform, in dein Programm. Der Rest ist einfach:

zeroes1 = [
    0.98811787286546693000, # comes from WinFilter
    0.00000000000000263287,
    2.96435361859640080000,
    0.00000000000000526575,
    2.96435361859640080000,
    0.00000000000000263287,
    0.98811787286546693000
]
poles1 = [
    1.00000000000000000000, # comes from WinFilter
    0.00000000013946281840,
    2.98743365009436030000,
    0.00000000021746471926,
    2.97494613251504750000,
    0.00000000001021009398,
    0.98751223597926108000
]

bandStopFromWinFilter = b2c_sft.Filter(zeroes1, poles1)
data = bandStopFromWinFilter.filterIterative(data)

Du siehst, dass du den Sinus bei 500 Hz auch mit einer Bandsperre aus WinFilter ausblenden kannst (leider nicht so exakt wie mit dem Notch. Der Notch ist Teil der Bibliothek, weil WinFilter diesen nicht anbietet, ich ihn aber für sehr sehr nützlich halte)

Mit diesen paar kleinen Hilfsmitteln kannst du schon die meisten Filteraufgaben lösen.

Und auf den ESP damit 🙂

Hier einfach ein Cheap-Yellow-Display als Referenzplatform

Benchmark

Als erstes brauchen wir ein Gefühl dafür wie lange die Berechnung einer FFT dauert. Dazu ein kleines Testprogramm:

import b2c_sft
import utime

lens = [16, 32, 64, 128, 256, 512, 1024, 2048]

for fftLen in lens:
    data = b2c_sft.generateSinus(fftLen, 100, 250, 2000)
    start = utime.ticks_us()
    b2c_sft.fft(data)
    time = utime.ticks_us() - start
    print("FFT-Länge " + str(fftLen) + " => " + str(time) + " us = " + str(time / 1000) + " ms") 

""" Ausgabe:
FFT-Länge 16 => 3476 us = 3.476 ms
FFT-Länge 32 => 5576 us = 5.576 ms
FFT-Länge 64 => 10649 us = 10.649 ms
FFT-Länge 128 => 21342 us = 21.342 ms
FFT-Länge 256 => 48962 us = 48.962 ms
FFT-Länge 512 => 97913 us = 97.91301 ms
FFT-Länge 1024 => 219556 us = 219.556 ms
FFT-Länge 2048 => 494714 us = 494.714 ms
"""

Sehr schön zu sehen ist, dass wir bei jeder Verdoppelung der FFT-Größe etwa doppelt so lange brauchen. Eine 512er FFT in ca. 100 ms find ich schon ganz anständig – damit kann man was anfangen 🙂

ESP-Programm

Das erste ESP-Programm hat nur etwa 50 Zeilen und zeigt die FFT mit dem lvgl-Chart an. Ist eine 256er FFT mit 4kHz Abtastrate. Man sieht, dass das ziemlich flüssig geht.

Hab ein Mikro an GPIO35 angeschlossen: https://www.amazon.de/dp/B08T1V3752 und etwas gepfiffen

#####################################################################
# Testprogramm für b2c_sft : box2code - Signal Filter Toolbox
# Stefan am 2024-01-05
#####################################################################
import machine, time
import uasyncio
import b2c_sft
import lvgl as lv
import display_driver

p0 = machine.Pin(21, machine.Pin.OUT)
p0 = 1

adc = machine.ADC(machine.Pin(35))
adc.atten(machine.ADC.ATTN_6DB)
adc.width(machine.ADC.WIDTH_9BIT)
tim = machine.Timer(1)
tsf = uasyncio.ThreadSafeFlag()

SAMPLE_FREQ = 4000
FFT_LEN = 256

data = []
timerStop = False

mid = 0
for i in range(0,1000):
    mid += adc.read()
mid = mid / 1000

def newAdcVal(p):
    global timerStop, tsf
    if timerStop:
        return
    data.append(200 * (adc.read()-mid))
    
    if len(data) == FFT_LEN:
        timerStop = True
        tsf.set()

chart = lv.chart(lv.scr_act())
chart.set_size(300, 220)
chart.align(lv.ALIGN.CENTER, 0, 0)
chart.set_range(lv.chart.AXIS.PRIMARY_Y, 0, 5000)
chart.set_style_size(0, lv.PART.INDICATOR)
ser = chart.add_series(lv.palette_main(lv.PALETTE.RED), lv.chart.AXIS.PRIMARY_Y)
chartData = []
for i in range(0, FFT_LEN // 2):
    chartData.append(0)
pcnt = len(chartData)
chart.set_point_count(pcnt)
chart.set_ext_y_array(ser, chartData)

def run():
    global timerStop
    while True:
        await tsf.wait()
        b2c_sft.hammingWindow(data)
        b2c_sft.fft(data)
        for i in range(0, len(data)):
            chartData[i] = int(data[i])
        chart.set_ext_y_array(ser, chartData)
        chart.refresh()
        data.clear()
        print("DONE")
        timerStop = False

tim.init(freq=SAMPLE_FREQ, mode=machine.Timer.PERIODIC, callback=newAdcVal)

uasyncio.run(run())

So wie filtert man denn dann jetzt ?

Dazu gehört jetzt nicht mehr viel. Definiere dir im ESP-Programm deinen Wunschfilter (hier wieder ein Notch) – Die Werte vom ADC müssen jetzt nur noch durch den Filter „durch“ – geschafft 🙂

notch = b2c_sft.Filter(b2c_sft.notchZeroes(500,SAMPLE_FREQ), b2c_sft.notchPoles(500,SAMPLE_FREQ))

def newAdcVal(p):
    global timerStop, tsf
    if timerStop:
        return
    data.append(notch.filterReal(adc.read() - 256))
    if len(data) == FFT_LEN:
        timerStop = True
        tsf.set()

Zum Ausprobieren einen Sinus von 0 bis 1000 Hz wobbeln – bei ca. 500 Hz sollte die Spitze verschwinden

Bibliothek

Die Bibliothek enthält natürlich wieder Bestandteile aus anderen Projekten:

• Hauptteil der FFT aus SimpleFFT – eine schöne kleine C++ Bibliothek, die ich schon in ein paar C++ Projekten benutzt habe
• Notch-Berechnung von hier
• Filter aus generiertem Code von WinFilter abgeleitet

#####################################################################
# b2c_sft : box2code.de - Signal Filter Toolbox
# Stefan am 2024-01-05
#####################################################################
import math

def isPowerOfTwo(num):
    return num and (not(num & (num - 1))) 

def rearrangeData(data):
    num_elements = len(data)
    buf = complex(0,0)
    target_index = 0
    bit_mask = 0
    for i in range(0, num_elements):
        if target_index > i:
            buf = data[target_index]
            data[target_index] = data[i]
            data[i]= buf
        bit_mask = num_elements
        while True:
            bit_mask = bit_mask >> 1
            if not (target_index & bit_mask):
                break
            target_index &= ~bit_mask
        target_index |= bit_mask

def makeTransform(data):
    counter = 0
    num_elements = len(data)
    local_pi = -math.pi
    nextbit = 0
    match = 0
    sine = 0.0
    delta = 0.0
    mult = complex(0,0)
    factor = complex(0,0)
    product = complex(0,0)

    i = 1
    while i < num_elements:
        nextbit = i << 1 
        delta = local_pi / i
        sine = math.sin(0.5 * delta)
        mult = complex(-2.0 * sine * sine, math.sin(delta))
        factor = complex(1.0)   

        for j in range(0,i):
            k = j
            while k < num_elements:
                counter += 1
                match = k + i
                product = data[match] * factor
                data[match] = data[k] - product
                data[k] += product
                k += nextbit
            factor = mult * factor + factor;
        i = i << 1
    return True
        
def backToReal(data):
    num_elements = len(data)
    for i in range(0, num_elements): 
        data[i] = abs(data[i]) / num_elements  
        
def fft(data):
    if not isPowerOfTwo(len(data)):
        print ("Not power of 2 !!!")
        return False
    rearrangeData(data)
    if not makeTransform(data):
        return False
    backToReal(data)
    n = len(data)
    for i in range(0, n // 2):
        data.pop()
    return True

def frequencyAxis(data, fs):
    num_elements = len(data)
    axis = []
    for i in range(0, num_elements):
        axis.append(fs / num_elements / 2 * i)
    return axis
    

def generateSinus(num_elements, amplitude, f, fs):
    sinus = []
    ta = 1 / fs
    t = 1/f
    for i in range(0, num_elements):
        sinus.append(amplitude * math.sin(i * 2*math.pi*ta/t))
    return sinus
    
def mixSignals(sig1, sig2):
    num_elements = min(len(sig1),len(sig2))
    data = []
    for i in range(0, num_elements):
        data.append(sig1[i] + sig2[i])
    return data


def notchPoles(f, fs):
    poles = []
    return poles
    
def notchZeroes(f, fs):
    zeroes = [0,0,0]
    zeroes[0] = 1.0
    zeroes[1] = -2.0 * math.cos(2 * math.pi * f / fs)
    zeroes[2] = 1.0
    return zeroes

def iirNotchZeroes(f, fs):
    b0=1.0
    b1=-2.0*math.cos(2.0*math.pi*f/fs);
    b2=1.0;
    return [b0,b1,b2]
    
def iirNotchPoles(f, fs):
    r=0.9;
    a0=1.0
    a1=-2.0*r*math.cos(2.0*math.pi*f/fs)
    a2=r*r
    return [a0,a1,a2]

class Filter:
    def __init__( self, zeroes, poles ):
        self.zeroes = zeroes
        self.poles = poles
        self.filterlen = len(zeroes)
        self.x = []
        self.y = []
        for i in range(0, self.filterlen + 1):
            self.y.append(0.0)
            self.x.append(0.0)  

    def filterReal(self, newSample): 
        # shift old samples
        for i in range(self.filterlen, 0, -1):
            self.x[i] = self.x[i-1]
            self.y[i] = self.y[i-1]
        # calc new output
        self.x[0] = newSample
        self.y[0] = self.zeroes[0] * self.x[0];
        for n in range(1, self.filterlen):
            self.y[0] += self.zeroes[n] * self.x[n]
            if len(self.poles) > 0:
                self.y[0] -= self.poles[n] * self.y[n]
        return self.y[0]
    
    def filterIterative(self, data):
        output = []
        for newSample in data:
            output.append(self.filterReal(newSample))    
        return output

def hammingWindow(data):
    for n in range(0, len(data)):
        data[n] = data[n] * (0.54 - 0.46 * math.cos((2 * math.pi / len(data)) * n))