FFT in Python kalibrieren/normieren

[Gustl B. (-gb-)]

Hallo,

ich verwende Python um Messdaten von ADCs zu erfassen. Aus diesen Daten 
möchte ich dann eine FFT rechnen lassen. Dasfunktioniert auch fein, ich 
sehe ein Spektrum im Plot, aber die Amplitude kann nicht stimmen.

Jetzt möchte ich das irgendwie kalibrieren oder normieren. Aber was 
sollte man verwenden? dBFS dBV oder was ganz anderes? Ich kann hier 
einen Sinus mit bekannter Amplitude anlegen und dann an den ADC Werten 
sehen wieviele Quantisierungsstufen verwendet wurden. Also ich weiß 
wieviele Volt einer Stufe entsprechen. Ich weiß aber nicht wie ich das 
der FFT mitteile.

Weil ich lieber Elektronik bastel als mich mit Software rumzuschlagen 
würde ich mich auch sehr über einen fertigen Codeschnipsel freuen.

Danke!

[void (Gast)]

Gustl B. schrieb:
> Ich weiß aber nicht wie ich das
> der FFT mitteile.

https://www.cbcity.de/die-fft-mit-python-einfach-erklaert

-> y-Axis: The Amplitude of the FFT Signal

[Analogopa (Gast)]

Gustl B. schrieb:
> Python um Messdaten von ADCs zu erfassen.
Da muss aber irgendwo noch eine HW dazwischen hängen. Wenn die nicht 
wichtig ist, dann ist der Umstand, dass es ADC-Werte sind auch nicht.

Bei einer Normierung hingegen wäre beides wichtig.

>Weil ich lieber Elektronik bastel als mich mit Software rumzuschlagen
>würde ich mich auch sehr über einen fertigen Codeschnipsel freuen.

Es wäre aber zweckmäßig, wenn du dich näher mit Signalverarbeitung und 
Messtechnik befassen würdest, um die Hintergründe zu verstehen.

>dBFS
Das ist so der einzige Weg.

Kalibrierung erfordert die gesamte Strecke.

[Gustl B. (-gb-)]

Analogopa schrieb:
> Da muss aber irgendwo noch eine HW dazwischen hängen.

Welche Hardware meinst Du denn genau? Der ADC ist Hardware, dahinter ist 
ein FPGA, der ist auch Hardware, dann kommt der FTDI USB Stein, auch 
Hardware, dann ein PC auch Hardware?! Aber die ändern an den ADC Werten 
nichts mehr.

Vor dem ADC hängt auch Hardware. Was muss ich denn davon wissen für die 
FFT?

Ich weiß folgendes:
1. Welche Spannung die von aussen angelegt wird welchem ADC-Wert 
entspricht.
2. Welche Spannung einer Quantisierungsstufe entspricht. Sprich die 
volle ADC Spanne entspricht x Volt, dann entspricht eine ADC Stufe x 
Volt geteilt durch die 2^n Stufen.
3. Kann ich einen externen Sinus anlegen mit bekannter Amplitude und 
dann eben an den ADC Werten sehen wie viele Stufen dieser Sinus 
verwendet.
4. Die Samplerate ist bekannt.

void schrieb:
> y-Axis: The Amplitude of the FFT Signal

Ja, wie das grundsätzlich funktioniert ist mir schon klar. Aber auch in 
deinem Link findet keine Normierung statt. Ausserdem ist da die 
Amplitude nicht in dB angegeben.

Analogopa schrieb:
>>dBFS
> Das ist so der einzige Weg.

Ist das wirklich so? Wenn ich weiß welche Spannung dem vollen ADC 
Wertebereich entspricht müsste doch auch dBV machbar sein?

Aber auch hier kommen Verwirrungen rein:
Der ADC liefert 16 Bit Werte als Offset Binary. Sprich von 0 bis 2^16-1. 
Setze ich dann als Referenz für dBFS 2^16 oder die Hälfte? Ich habe mal 
Beispielcode gesehen der für ein 16 Bit Signal das aber symmetrisch um 
die 0 ist nur 2^15 für die Referenz verwendet hat.

[-gb- (Gast)]

So, habe mal ein kleines Beispiel gebaut, den Großteil habe ich aus im 
Internet, aber es funktioniert und erzeugt einen Plot. Im Anhang ist 
jetzt das Python und eine Datei die Samplewerte enthält. Die kommen von 
einem 16 Bit ADC und wurden mit 25 MSamples/s aufgenommen.

import numpy
import matplotlib.pyplot as plt
from scipy.fftpack import fft, fftshift
from scipy.signal import blackman

ADC_FS = 25.0 #MHz
ADC_Bits = 16
ADC_Steps = 2**ADC_Bits
ADC_N_Samples = 2**13
ADC_Samples = []
file = open("sine64k_10MHz.txt", "r")
# Read Samples
linecounter = 0
for line in file:
  if linecounter < ADC_N_Samples:
    ADC_Samples.append(int(line))
    linecounter +=1
file.close()

freq = numpy.linspace(-ADC_FS/2, ADC_FS/2, ADC_N_Samples)
A = fft(ADC_Samples*blackman(ADC_N_Samples))
response = 20 * numpy.log10(numpy.abs(fftshift(A / abs(A[1:ADC_N_Samples]).max())))
plt.plot(freq, response)
plt.axis([-0.5, 0.5, -140, 0])
plt.xlim(xmax = ADC_FS/2, xmin = 0)
plt.title("Blackman window")
plt.ylabel("Amplitude [dB]")
plt.xlabel("f [MHz]")
plt.show()

[-gb- (Gast)]

Ich schrieb es funktioniert, aber damit meinte ich nicht, dass es die 
Lösung ist. Die Amplitude stimmt weiterhin nicht.

[void (Gast)]

Gustl B. schrieb:
> Ja, wie das grundsätzlich funktioniert ist mir schon klar. Aber auch in
> deinem Link findet keine Normierung statt.

Doch genau das. Für die x- (Frequenz) und die y- Achse (Amplitude). Mit 
schrittweise erklärtem Python Code. Hätte man aber lesen und vertehen 
müssen. Einfach nur den Großteil aus dem Internet abschreiben 
funktioniert halt nicht.

Gustl B. schrieb:
> Ausserdem ist da die
> Amplitude nicht in dB angegeben.

Dann multipliziert man sie noch einmal mit 20*log(n) und ist fertig.

[Gustl B. (-gb-)]

Nun, auf der hier empfohlenen Seite wird aus einem Sinus von -100 bis 
+100, also mit Amplitude 200 ein Peak im Spektrum mit Amplitude von 200.

Ich habe bisher Sampledaten von einem Fullscale Sinus erzeugt und 
Rauschen draufaddiert. Der Sinus hat 5 MHz und die Abtastrate beträgt 25 
MHz. Die Samples haben Werte von 0 bis 2**16-1. Das zeigt das Bild 
FFT_werte.png.

Dann habe ich das mit einer Fensterfunktion multipliziert.

Y = np.hanning(len(werte))*werte

Das sieht jetzt schon anders aus als auf der Webseite, weil mein Signal 
nicht von -2**15 nach +2**15 geht sondern nur positiv ist. Aber das ist 
eben dann auch das was ich später vom AD Wandler bekomme, da geht auch 
DC durch und das kann an beliebiger Stelle liegen. Raus kommt 
FFT_werte_fenster.png.

Das wurde dann in die FFT gesteckt

Y = np.fft.fft(Y)
Y = 2.0*np.abs(Y[:N])/N

und es entsteht ein Peak. Der hat aber nur die halbe Höhe vom vollen 
Wertebereich. Warum weiß ich nicht. Zu sehen in FFT_peak.png.

Dann wurde mir noch folgendes empfohlen:
void schrieb:
> Dann multipliziert man sie noch einmal mit 20*log(n) und ist fertig.
Unklar bleibt was n hier ist. Die Anzahl der Samples? Mein Y?

Verwende ich

Y = 20*np.log(Y)

 kommt etwas sinnloses bei raus. Zumindest machen die Werte keinen Sinn. 
FFT_peak_dB.png

Wie geht das, dass ein Peak mit maximal möglicher Höhe bei 0 dBFs liegt?

[Achim S. (Gast)]

Gustl B. schrieb:
> Unklar bleibt was n hier ist. Die Anzahl der Samples? Mein Y?

Dein Y.

Gustl B. schrieb:
> Verwende ichY = 20*np.log(Y) kommt etwas sinnloses bei raus. Zumindest
> machen die Werte keinen Sinn.

du sollst auch den Zehner-Logarithmus nehmen, nicht den natürlichen 
Logarithmus.

https://docs.scipy.org/doc/numpy/reference/generated/numpy.log.html

Gustl B. schrieb:
> Wie geht das, dass ein Peak mit maximal möglicher Höhe bei 0 dBFs liegt?

rein phänomenologisch: dividiere dein Y vor der Logarithmierung durch 
den Wert, den du als maximalen Peak rausgemessen hast (also rund 2^15).

Wenn es um Verständnis geht: schau dir die Dokumentation deiner Funktion 
an
https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fft.html
und die zugehörigen implementation details
https://docs.scipy.org/doc/numpy/reference/routines.fft.html#module-numpy.fft

[Gustl B. (-gb-)]

Achim S. schrieb:
> du sollst auch den Zehner-Logarithmus nehmen, nicht den natürlichen
> Logarithmus.

Ok, da kommt dann FFT_peak_dB_log10.png bei raus.

Achim S. schrieb:
> rein phänomenologisch: dividiere dein Y vor der Logarithmierung durch
> den Wert, den du als maximalen Peak rausgemessen hast (also rund 2^15).

Fein, da kommt dann was bei rum was für mich ok aussieht. Danke!

Edit: Jetzt habe ich ein zweites Signal draufaddiert bei 10 MHz mit 
1/1000 der Amplitude. Und siehe da, es liegt wunderbar bei -60 dB. Ich 
bin sehr zufrieden, vielen Dank!

So, und zum Schluss noch Plots von echten Messwerten. Die 10 MHz 
Referenz ist ein TSW2110EVM.

[Tobias (. (Gast)]

Achim S. schrieb:
> rein phänomenologisch: dividiere dein Y vor der Logarithmierung durch
> den Wert, den du als maximalen Peak rausgemessen hast

Das ist aber nicht perfekt. Der ermittelte Peak varriert doch mit der 
Abdeckung der FFT-Frequenz / der Breite.

?

[Gustl B. (-gb-)]

Also ich habe durch 2^15 dividiert. Ich verstehe aber nicht wieso. Der 
ADC löst ja in 2^16 Stufen auf. Kommt der Faktor 1/2 daher, weil man von 
der FFT auch nur eine Hälfte verwendet?

In der nächsten Revision werde ich statt des einen ADA4932-2 zwei 
ADA4945 verwenden und dann zwischen diesen und dem ADC noch einen 
HMC960LP4E verbauen. Den habe ich heute auf einer Testplatine über SPI 
bespaßt und ich bin sehr zufrieden. Klar wird er Rauschen hinzufügen, 
aber das ist mir bei diesem Lernbastelprojekt egal.

[Lukas (Gast)]

Gustl B. schrieb:
> Also ich habe durch 2^15 dividiert. Ich verstehe aber nicht wieso. Der
> ADC löst ja in 2^16 Stufen auf. Kommt der Faktor 1/2 daher, weil man von
> der FFT auch nur eine Hälfte verwendet?

Wenn du einen Sinus mit 16 Bit auflöst, wie groß kann dann die max. 
Amplitude(!) sein?

[Gustl B. (-gb-)]

Stimmt, daran hatte ich nicht gedacht.