Re: protonovy magnetometr, slaby signal
Vojtěch Petrucha
petrucha na volny.cz
Neděle Květen 30 12:41:05 CEST 2021
já do tohohle moc nešťoural, tak mohu jen zkopírovat co má k tomu student v dp a přiložit příslušný python kód.
našel jsem ty fft spektra v příloze dp, přidal jsem na https://maglab.fel.cvut.cz/proton-mag-dalsi-pokusy/ <https://maglab.fel.cvut.cz/proton-mag-dalsi-pokusy/>
asi by šlo použít jen 400k vzorků místo 4M a nějak to interpolovat.. asi vyzkoušíme
v.
K analýze měřených dat se také velmi často využívá analytického signálu, který je reprezentován
pomocí komplexních čísel. Reálnou část tvoří skutečná naměřená data a imaginární část Hilbertova
transformace naměřených dat. Matematicky lze takový signál vyjádřit pomocí vztahu:
𝑠𝑎 (𝑡) = 𝐴(𝑡)𝑒𝑗𝜑(𝑡) = 𝐴(𝑡) ⋅ cos(𝜑(𝑡)) + 𝐴(𝑡) ⋅ 𝑗 sin(𝜑(𝑡)). (2.6)
Pro bližší představu můžeme na obrázku 2.3 vidět grafické znázornění tří period harmonického signálu
převedeného na analytický signál. Zkonstruvaný signál si tedy lze představit jako pohyb popisující
jednotkovou kružnici v čase. Díky tomu můžeme snadno určit okamžitou amplitudu a frekvenci, a právě
proto se této techniky využívá při frekvenční analýze [7]. V našem případě je možné tuto metodu využít
k určení průměrné frekvence a následně magnetické indukce odpovídající měřenému precesnímu
signálu.
[7] B. Boashash, „Estimating and interpreting the instantaneous frequency of a signal. I. Fundamentals,“
Proceedings of the IEEE, sv. 80, č. 4, pp. 520 - 538, 1992 .
def computeFrequencyWithHilbertTransform(signal, startSample, endSample, fs):
selectedSignal = []
d_angle = []
d_freq = []
skippedEnd = 50
hilbertSignalLength = endSample - startSample - skippedEnd
n = endSample - startSample - skippedEnd - 1
for k in range(len(signal)):
if((k >= startSample) and (k < endSample)):
selectedSignal.append(signal[k] - 32768) #shift signal to 0
analytic_signal = hilbert(selectedSignal)
#omit last 50samples of analytic signal
selectedSignal = selectedSignal[:len(selectedSignal)-skippedEnd]
analytic_signal = analytic_signal[:len(analytic_signal)-skippedEnd]
#d_angle = np.unwrap(np.angle(analytic_signal))
#determine instantenous angle
skippedBegining = 20
for k in range(len(analytic_signal)-1):
if(k > skippedBegining):
d_angle.append(np.angle(analytic_signal[k+1]/analytic_signal[k]))
#determine instantenous frequency
for k in range(len(d_angle)):
d_freq.append(d_angle[k] / (1/fs)/(2*math.pi))
d_angle = np.cumsum(d_angle)
# dtermine average frequency
diff_d_angle = max(d_angle) - min(d_angle)
tau=2 * math.pi * ((n-skippedBegining)/fs) / diff_d_angle
f_average = 1/tau
#plot - hilbert transform
if(enablePlot_hilbert):
fig = plt.figure(figsize=(30, 12))
plt.subplot(2, 1, 1)
plt.plot(signal, label='Signál')
plt.plot(np.imag(analytic_signal), label='Hilbertova transformace naměřených dat')
plt.xlabel("Vzorky")
plt.ylabel("Amplituda")
plt.legend()
plt.title(graphTitle + " Hilbertova transformace")
plt.subplot(2, 4, 5)
plt.xlim(0, 400)
plt.plot(signal)
plt.plot(np.imag(analytic_signal))
plt.xlabel('step')
plt.ylabel('Amplitude')
plt.subplot(2, 4, 6)
plt.xlim(hilbertSignalLength/3 - 200, hilbertSignalLength/3 + 200)
plt.plot(signal)
plt.plot(np.imag(analytic_signal))
plt.xlabel('step')
plt.ylabel('Amplitude')
plt.subplot(2, 4, 7)
plt.xlim(hilbertSignalLength/3 * 2 - 200, hilbertSignalLength/3 * 2 + 200)
plt.plot(signal)
plt.plot(np.imag(analytic_signal))
plt.xlabel('step')
plt.ylabel('Amplitude')
plt.subplot(2, 4, 8)
plt.xlim(hilbertSignalLength - 400, hilbertSignalLength)
plt.plot(signal)
plt.plot(np.imag(analytic_signal))
plt.xlabel('step')
plt.ylabel('Amplitude')
fig.savefig('generatedFiles/figures/' + graphTitle + '_hilbert.png', dpi=fig.dpi)
plt.close(fig)
# print(" measured signal and hilbert transform")
# plot reconstructed signal vs measured signal
if(enablePlot_reconstructed):
#create reconstructed signal with detemined frequency
reconstructedSig = np.zeros(endSample- startSample)
amplitude =max(np.abs(signal)*0.5)
#omit also last 50 samples
for k in range(endSample- startSample):
reconstructedSig[k] = (amplitude*math.cos(2*math.pi*f_average*(k)/fs))
#plot
fig = plt.figure(figsize=(30, 12))
plt.subplot(2, 1, 1)
plt.plot(selectedSignal, label='Naměřený signál')
plt.plot(reconstructedSig, label='Zrekonstruovaný signál')
plt.xlabel("Vzorky")
plt.ylabel("Amplituda")
plt.legend()
plt.title(graphTitle + " Rekonstrukce signálu")
plt.subplot(2, 4, 5)
plt.xlim(0, 400)
plt.plot(selectedSignal)
plt.plot(reconstructedSig)
plt.xlabel('step')
plt.ylabel('Amplitude')
plt.subplot(2, 4, 6)
plt.xlim(hilbertSignalLength/3 - 200, hilbertSignalLength/3 + 200)
plt.plot(selectedSignal)
plt.plot(reconstructedSig)
plt.xlabel('step')
plt.ylabel('Amplitude')
plt.subplot(2, 4, 7)
plt.xlim(hilbertSignalLength/3 * 2 - 200, hilbertSignalLength/3 * 2 + 200)
plt.plot(selectedSignal)
plt.plot(reconstructedSig)
plt.xlabel('step')
plt.ylabel('Amplitude')
plt.subplot(2, 4, 8)
plt.xlim(hilbertSignalLength - 400, hilbertSignalLength)
plt.plot(selectedSignal)
plt.plot(reconstructedSig)
plt.xlabel('step')
plt.ylabel('Amplitude')
fig.savefig('generatedFiles/figures/' + graphTitle + '_reconstructed.png', dpi=fig.dpi)
plt.close(fig)
# print(" measured signal and reconstructed signal")
# plot instantaneous angle
if(enablePlot_instantaneousAngle):
fig = plt.figure(figsize=(30, 12))
plt.subplot(2, 1, 1)
plt.plot(d_angle)
plt.xlabel("Vzorky")
plt.ylabel("Úhel [rad]")
plt.title(graphTitle + " Okamžitý úhel")
plt.subplot(2, 4, 5)
plt.xlim(0, 400)
plt.plot(d_angle)
plt.xlabel("Vzorky")
plt.ylabel("Úhel [rad]")
plt.subplot(2, 4, 6)
plt.xlim(hilbertSignalLength/3 - 200, hilbertSignalLength/3 + 200)
plt.plot(d_angle)
plt.xlabel("Vzorky")
plt.ylabel("Úhel [rad]")
plt.subplot(2, 4, 7)
plt.xlim(hilbertSignalLength/3 * 2 - 200, hilbertSignalLength/3 * 2 + 200)
plt.plot(d_angle)
plt.xlabel("Vzorky")
plt.ylabel("Úhel [rad]")
plt.subplot(2, 4, 8)
plt.xlim(hilbertSignalLength - 400, hilbertSignalLength)
plt.plot(d_angle)
plt.xlabel("Vzorky")
plt.ylabel("Úhel [rad]")
fig.savefig('generatedFiles/figures/' + graphTitle + '_instantaneousAngle.png', dpi=fig.dpi)
plt.close(fig)
# print(" instantaneous angle")
# plot instantaneous frequency
if(enablePlot_instantaneousFrequency):
fig = plt.figure(figsize=(30, 12))
plt.subplot(2, 1, 1)
plt.plot(d_freq)
plt.ylim(1800, 2400)
plt.xlabel("Vzorky")
plt.ylabel("Frekvence [Hz]")
plt.subplot(2, 4, 5)
plt.xlim(0, 400)
plt.ylim(1800, 2400)
plt.plot(d_freq)
plt.xlabel("Vzorky")
plt.ylabel("Frekvence [Hz]")
plt.subplot(2, 4, 6)
plt.xlim(hilbertSignalLength/3 - 200, hilbertSignalLength/3 + 200)
plt.ylim(1800, 2400)
plt.plot(d_freq)
plt.xlabel("Vzorky")
plt.ylabel("Frekvence [Hz]")
plt.subplot(2, 4, 7)
plt.xlim(hilbertSignalLength/3 * 2 - 200, hilbertSignalLength/3 * 2 + 200)
plt.ylim(1800, 2400)
plt.plot(d_freq)
plt.xlabel("Vzorky")
plt.ylabel("Frekvence [Hz]")
plt.subplot(2, 4, 8)
plt.xlim(hilbertSignalLength - 400, hilbertSignalLength)
plt.ylim(1800, 2400)
plt.plot(d_freq)
plt.xlabel("Vzorky")
plt.ylabel("Frekvence [Hz]")
fig.savefig('generatedFiles/figures/' + graphTitle + '_instantaneousFreq.png', dpi=fig.dpi)
plt.close(fig)
return f_average
______________________________________________________________
> Od: "balu" <balu na k-net.fr>
> Komu: hw-list na list.hw.cz
> Datum: 30.05.2021 12:05
> Předmět: Re: protonovy magnetometr, slaby signal
>
Este jedna otazka k tomu Hilbertovi. Da sa strucne popisat ta metoda, aby som si urobil predstavu?
Podobny typ problemu riesim na 1 Gb/s signaloch z LHC, mozno by sa dalo nieco aplikovat aj na toto.
b.
On 30/05/2021 11:53, Vojtěch Petrucha wrote:zdravim,
ideálně se hledá jeden bin s maximální hodnotou, přičemž přibližná frekvenční hodnota toho binu může být známá (z měření frekvence komparátorem..), možná by šlo použít i data třeba z okolních binů toho maximálního a ještě to interpolovat..
student používal na pc skript v python a knihovnu scipy.fft.. a předpokládám že ve výsledku prostě hledal maximum, nějak to tam teď nevidím, musím se zeptat..
v.
______________________________________________________________
> Od: "balu" <balu na k-net.fr> <balu na k-net.fr>
> Komu: hw-list na list.hw.cz <hw-list na list.hw.cz>
> Datum: 29.05.2021 17:12
> Předmět: Re: protonovy magnetometr, slaby signal
>
v spektre sa hlada jedna dominantna ciara, ktorej frekvenciu priblizne pozname? Alebo je to vseobecne spektrum, ktore moze obsahovat hocico?
b.
On 29/05/2021 14:20, Vojtěch Petrucha wrote:zdravim,
tak jsme našli ještě trochu času a asi už to funguje, viz poslední obrázky na https://maglab.fel.cvut.cz/proton-mag-dalsi-pokusy/ <https://maglab.fel.cvut.cz/proton-mag-dalsi-pokusy/>
změřili jsme zánik proudu pro různé hodnoty toho tlumícího rezistoru, pro 1k to vypadalo velmi rozumně, tak jsme zkrátili dobu kdy je sepnutý, aby více zbytečně netlumil, vyběhli před FEL na trávník a signál prakticky totožný s tím s přepináním pomocí relé...
takže moc díky za pomoc
chtěl bych zkusit to zjišťování frekvence pomocí fft, výsledky vypadaly zajímavě. k dispozici 30-40k vzorků při 44.1kSa/s , potřeba alespoň 0.01 Hz rozlišení, to vede na 4M (2^22) délku fft, ale předpokládám že výpočet lze nějak zjednodušit, když by 99% byly nuly... ideální by bylo počítat přímo v mcu (stm32f767vi, 512kByte ram...) v době následující polarizační fáze - 3s k dispozici... ale i připlácnutí raspberry pi zero by asi také šlo...
v.
----------
_______________________________________________
HW-list mailing list - sponsored by www.HW.cz <http://www.HW.cz>
Hw-list na list.hw.cz <Hw-list na list.hw.cz>
http://list.hw.cz/mailman/listinfo/hw-list <http://list.hw.cz/mailman/listinfo/hw-list>
_______________________________________________HW-list mailing list - sponsored by www.HW.cz <http://www.HW.cz>Hw-list na list.hw.cz <Hw-list na list.hw.cz>http://list.hw.cz/mailman/listinfo/hw-list <http://list.hw.cz/mailman/listinfo/hw-list>
----------
_______________________________________________
HW-list mailing list - sponsored by www.HW.cz
Hw-list na list.hw.cz
http://list.hw.cz/mailman/listinfo/hw-list <http://list.hw.cz/mailman/listinfo/hw-list>
------------- další část ---------------
HTML příloha byla odstraněna...
URL: <http://list.hw.cz/pipermail/hw-list/attachments/20210530/50c971e5/attachment.html>
Další informace o konferenci Hw-list