10 Hz in het oor

[Xilvo]

Dit krijg je (f₀=10 Hz, samengesteld t.m. de duizendste harmonische, links zonder, rechte met willekeurige fase.
blok10Hz.png

De vorm blijft hier wel redelijk goed behouden.

Meer dan redelijk, sterker, die is perfect periodiek.
Dat is logisch, ik transformeer niet maar genereer het signaal direct door de verschillende harmonische bij elkaar op te tellen. De faseverschuiving ligt per harmonische vast, dan moet het signaal wel periodiek worden.

[wnvl1]

1024, opnieuw niet perfect, maar wel redelijk.
Ook veel meer vervormd dan jullie figuren.
De reconstructie werkt perfect als ik niets verander aan de fase.
De fase veranderen doe ik door te vermenigvuldigen met een willekeurig punt op de eenheidscirkel in het complexe vlak voor de eerste helft en dan te vermenigvuldigen met het toegevoegde punt in de tweede helft.

Code: Selecteer alles

import numpy as np
import scipy 
from scipy import signal
import matplotlib.pyplot as plt

n=1024
t = np.linspace(0, 1, n)

zaagtand = signal.square(2 * np.pi * 5 * t)
plt.plot(t, zaagtand)
plt.title("blokgolf")

zaagtand_fft = scipy.fft.fft(signal.sawtooth(2 * np.pi * 5 * t)) 
random_fases = np.exp(1j *np.random.rand(n // 2 + 1) * 2 * np.pi)  
random_fases_gespiegeld = np.concatenate([random_fases, np.conj(random_fases[-2:0:-1])])
zaagtand_fft_gewijzigde_fase = zaagtand_fft * random_fases_gespiegeld
zaagtand_reconstructie = scipy.fft.ifft(zaagtand_fft_gewijzigde_fase)
plt.plot(t, zaagtand_reconstructie)
plt.title("blokgolf willekeurige fase")

[wnvl1]

Meer dan redelijk, sterker, die is perfect periodiek.
Dat is logisch, ik transformeer niet maar genereer het signaal direct door de verschillende harmonische bij elkaar op te tellen. De faseverschuiving ligt per harmonische vast, dan moet het signaal wel periodiek worden.

Mijn opmerking was anders bedoeld. Ik had het niet over de periodiciteit, maar over de gelijkenis in vorm tussen één blok en wat jij bekomt na transformatie van de fases.

[Xilvo]

Mijn opmerking was anders bedoeld. Ik had het niet over de periodiciteit, maar over de gelijkenis in vorm tussen één blok en wat jij bekomt na transformatie van de fases.

Dat is toeval, dat hangt voornamelijk af van welke faseverschuiving de lagere harmonischen (derde, vijfde) krijgen.
Een andere keer, met andere willekeurige faseverschuivingen, ziet het er zo uit:

[HansH]

Zelfde code, maar nu met een sampling van 1000 punten ipv 100punten.

bedoel je van 100 naar 1000 samples over de periode van 1s?
met 100 samples per 1 seconde is de sample afstand 10ms dus dan kun je frequentie componenten van 50 Hz of meer al niet meer weergeven dus als je die optelt krijg je fouten.
wat FFT daarmee te maken heeft zie ik nog even niet want je bent voor de signaal reconstructie sinussen aan het optellen. met FFT bereken je frequentiecomponenten, maar dat is hier niet nodig want je weet immers al uit de theorie welke componenten dat zijn en welke amplitude en fase.

[wnvl1]

bedoel je van 100 naar 1000 samples over de periode van 1s?

ja

met 100 samples per 1 seconde is de sample afstand 10ms dus dan kun je frequentie componenten van 50 Hz of meer al niet meer weergeven dus als je die optelt krijg je fouten.

maar de periodiciteit zou behouden moeten blijven

wat FFT daarmee te maken heeft zie ik nog even niet want je bent voor de signaal reconstructie sinussen aan het optellen. met FFT bereken je frequentiecomponenten, maar dat is hier niet nodig want je weet immers al uit de theorie welke componenten dat zijn en welke amplitude en fase.

FFT is een 'luie' oplossing. Als je FFt gebruikt hoef je zelf niet meer de integralen uit te rekenen om de Fourier reeks op te stellen. In scipy zijn alle klassieke functies uit de signaalverwerking ingebouwd. Overigens wat een FFT doet is in de grond de Fourierreeks uitrekenen.

[HansH]

maar de periodiciteit zou behouden moeten blijven

dat klopt. je zou eens kunnen kijken als je alleen 1 of enkele van de frequentie componenten bekijkt als deel van de sommatie. dan moet je kunnen zien waar het probleem ontstaat.

[HansH]

FFT is een 'luie' oplossing. Als je FFt gebruikt hoef je zelf niet meer de integralen uit te rekenen om de Fourier reeks op te stellen. In scipy zijn alle klassieke functies uit de signaalverwerking ingebouwd. Overigens wat een FFT doet is in de grond de Fourierreeks uitrekenen.

dat is duidelijk, maar ik bedoelde dat de uitkomt van die integralen al gedaan is voor jou en je zo de formule voor elke frequentiecomponent beschikbaar hebt van de standaard vormen zoals blok, zaag en driehoek. (zoals xilvo ook al aangaf) zo heb ik het ook in LTspice grstopt uitgaande van de formule voor een blok

[wnvl1]

Mijn insteek is dat ik graag dit soort van problemen oplos met als uitgangspunten (1) zo snel mogelijk/zo weinig mogelijk lijnen code, (2) zo flexibel/algemeen mogelijk, (3) gebruik maken van alle libraries die al beschikbaar zijn en (4) als ik leuke wiskunde kan toepassen, kan dat het vorige zeker overrulen.

Bij deze opgave gingen mijn gedachte daarom onmiddellijk uit naar FFT.

[HansH]

dus je gebruikt de FFT om de frequentiecomponenten te berekenen en daarna tel je die componenten als functie van de tijd (of sample momenten in jouw geval) op om weer het oorspronkelijke signaal te resonstrueren als ik het goed begrijp. en als alternatief dezelfde frequentie componenten die uit de FFT rollen, maar dan met een random fase voor elke component?

[wnvl1]

Dat klopt. Dit zijn de vijf lijnen code waarin dit gebeurt.

Code: Selecteer alles

zaagtand_fft = scipy.fft.fft(signal.sawtooth(2 * np.pi * 5 * t)) 
random_fases = np.exp(1j *np.random.rand(n // 2 + 1) * 2 * np.pi)  
random_fases_gespiegeld = np.concatenate([random_fases, np.conj(random_fases[-2:0:-1])])
zaagtand_fft_gewijzigde_fase = zaagtand_fft * random_fases_gespiegeld
zaagtand_reconstructie = scipy.fft.ifft(zaagtand_fft_gewijzigde_fase)

[Xilvo]

De IFFT geeft complexe waardes, die plot je vervolgens. Gaat dat wel goed?

[HansH]

De IFFT geeft complexe waardes, die plot je vervolgens. Gaat dat wel goed?

dat zijn als het goed is 2 toegevoegd complexe waardes die samen een reele amplitude en fase geven.
hier een voorbeeldje met fourier.

fourierheentrug

(455.25 KiB) 28 keer gedownload

[wnvl1]

Het complexe deel is heel klein, dat zou moeten vallen in de categorie numerieke fouten. Ik heb bvb na IFFT

[0.64439172-0.00197186j 0.76540039-0.0005845j 0.55170809-0.00197186j ...
0.69290493-0.0005845j 0.75350153-0.00197186j 0.4890989 -0.0005845j ]

Door van de random fase factoren in het tweede deel van de array die uit de FFT komt het complex toegevoegde te nemen, behoud je 'in principe' een reëel signaal. Dat gaat natuurlijk nooit helemaal exact kloppen.

[Xilvo]

Het complexe deel is heel klein, dat zou moeten vallen in de categorie numerieke fouten. Ik heb bvb na IFFT

[0.64439172-0.00197186j 0.76540039-0.0005845j 0.55170809-0.00197186j ...
0.69290493-0.0005845j 0.75350153-0.00197186j 0.4890989 -0.0005845j ]

Door van de random fase factoren in het tweede deel van de array die uit de FFT komt het complex toegevoegde te nemen, behoud je 'in principe' een reëel signaal. Dat gaat natuurlijk nooit helemaal exact kloppen.

Maar dit klopt wel erg slecht, dat zijn geen afrondingsfouten meer.

Component nul van de transform is het gemiddelde, daar is fase betekenisloos.
Component n/2 moet zijn eigen geconjugeerde zijn, moet reëel zijn.
Als ik beide componenten in random_fases_gespiegeld 1 maak, dan krijg ik in ieder geval een IFFT waarbij de imaginaire delen allemaal echt klein zijn, meestal kleiner dan 1E-15.