verauschtes Signal

[fourier (Gast)]

Hallo ich habe ein verauschtes, periodisches Signal und möchte das 
Rauschen dort herausfiltern mithilfe der Fouriertransformation und 
inverse fft.

Matlab Code:

fourier_samples = load('data_ha.mat');
four_sa= fourier_samples.data;

Fs = 4347828; %Abtastfrequenz

plot(four_sa)

[y_freq,freq_rng] = pos_fft(four_sa, Fs);
figure
plot(freq_rng,abs(y_freq))
xlim([0,10000])

[pks,loc] = findpeaks(abs(y_freq), 'MinPeakHeight', 0.001, 'MinPeakDistance', 0.15);


y_freq(1:241) = 0+0*1i;
y_freq(243:723) = 0+0*1i;
y_freq(725:1206) = 0+0*1i;
y_freq(1208:1688) = 0+0*1i;
y_freq(1670:2170) = 0+0*1i;
y_freq(2172:end-2172) = 0+0*1i;
y_freq(end-240:end) = 0+0*1i;
y_freq(end-722:end-242)=  0+0*1i;
y_freq(end-1205:end-724)=  0+0*1i;
y_freq(end-1687:end-1207)=  0+0*1i;
y_freq(end-2169:end-1689)=  0+0*1i;

y_freq_1 = reshape(y_freq,1,length(y_freq));
y_freq_2 = flip(y_freq);
y_freq_2 = reshape(y_freq_2,1,length(y_freq_2));

y_freq = [y_freq_1 y_freq_2];

figure


new = ifft(y_freq);
new = real(new);

plot(new)





Funktion:

function[X,freq] = pos_fft(x,Fs)

N = length(x);
k = 0:N-1;
T = N/Fs;
freq = k/T;
X = fft(x)/N;
cutOff = ceil(N/2);
X = X(1:cutOff);
freq = freq(1:cutOff);

Leider bekomme ich nicht das gewünschte Ergebnis, obwohl mit das 
Amplitudenspektrum ausgeben lasse und mir die einzelnen Peaks 
rausfiltere .
Um dann die restlichen Frequenzen auf Null zu setzen und ifft 
anzuwenden, damit ich anschließend das Signal ohne Rauschen darstellen 
kann.
Woran könnte das liegen?

MFG

[Markus B. (russenbaer)]

Servus,

Soweit ich das im Code erkennen kann und unter der Annahme das es sich 
beim Ursprungssignal um ein reelles Signal handelt (zumindest setzt Du 
das mit pos_fft voraus), setzt Du das Spektrum für die Synthese falsch 
zusammen.
Mir geht das Konjugieren ab...
Weiters glaube ich das Du einen Indexfehler hast (zumindest wenn die FFT 
geradzahlig ist [Der NQ bin kommt doppelt vor und den bin bei der 
Abtastfrequenz brauchst auch nicht.]).

Zum Code: bitte mach Kommentare...

lg
Russenbaer

[fourier (Gast)]

Danke für die Tipps.

Ich habe mich dazu entschieden einfach nochmal anzufangen und step by 
step voranzukommen.


folgenden Code habe ich bisher geschrieben:

fourier_samples = load('data_ha.mat');
four_sa= fourier_samples.data;


Fs = 4347828; %Abtastfrequenz

plot(four_sa) %Darstellung des verauschten Signals

Y = fft(four_sa); %Diskrete Fouriertransformation angewandt auf das Ursprungssignal

n=length(Y)/2; %Länge halbieren um nur die positive Seite des Spektrums darzustellen
fq = 0:n-1; % frequenzbins
fq = fq * Fs /length(four_sa); %für Darstellung in Hz
figure(2)
plot(fq(1:10000),abs(Y(1:10000))/length(four_sa)/2) %Darstellung des Betragsspektrum im  Frequenzbereich bis 10 KHz

Im Anhang liegt ein JPG mit dem Amplitudenspektrum bis 10 KHz.


Jetzt bestimmt man mithilfe der Grafik und der Variable 'Y' die Peaks 
und
kann ifft benutzen. Für diese Methode muss ich doch alle Komponenten die 
kein Peak sind Null setzen und dann ifft anwenden.
Müssen die gespiegelten Frequenzpeaks auch Null gesetzt werden ?
Anschließend nimmt man sich die Realteile und hat das gefilterte Signal?



Ein andere Möglichkeit wäre sich die Peaks zu nehmen und damit eine 
Summe aus Harmonischen(Sin/Cos) zu konstruieren. Nur da weiß ich nicht 
wie hoch die Frequenzen und Koeffizienten sind.

Sind die Koeffizienten die Absolutbeträge der Peaks?
Sind die Frequenzen für die Harmonischen die Faktoren von der Abszisse 
meines Plots von dem Amplitudenspektrum?

MFG

[Markus B. (russenbaer)]

Servus,

Für ein reelles Signal gilt X(-jw) = X*(jw), d.h. das Spektrum ist 
konjugiert komplex. Siehe z.B. 
https://www.nt.tuwien.ac.at/wp-content/uploads/2016/01/Doblinger_SuS_3Ed_print.pdf 
S.70

D.h. bei einer geraden FFT entspricht Nyquist-bin und der DC-bin einem 
Spiegelpunkt im konjugiert komplexen Sinn.
z.B. brauchst Du für eine FFT der Größe 8 eines reelwertigen Signals im 
Frequenzbereich nur 5 bins: bin 0 is bin 4
Diese kannst Du modifizieren. Für die Synthese musst Du für die IFFT 
diese bins wieder generieren.

B(5) = conj(B(3))
B(6) = conj(B(2))
B(7) = conj(B(1))

Damit hast Du für die IFFT alle bins die du für die Rückwandlung 
brauchst - bin 0 bis bin 7.

Je nach Definition des Paares FFT/IFFT brauchst du noch den Faktor 1/N 
entweder bei der FFT, IFFT oder 1/sqrt(N) bei beiden. Das musst Du Dir 
in der Matlab Hilfe ansehen.

"Müssen die gespiegelten Frequenzpeaks auch Null gesetzt werden ?
Anschließend nimmt man sich die Realteile und hat das gefilterte 
Signal?"
Nein und noch einmal Nein


lg
RB

[Sven B. (scummos)]

Das klappt so eher nicht. Du willst vermutlich einen FIR- oder 
IIR-Filter hoher Ordnung synthetisieren, der dann auch für das Signal 
einen sinnvollen Phasenverlauf hat.

Was willst du denn erreichen? Das Spektrum sieht jedenfalls schon 
ziemlich sauber aus.

[fourier (Gast)]

Mein Ziel ist es das Rauschen aus dem Signal zu bekommen um das 
periodische Signal zu bekommen.
Man sieht an meiner Abbildung ja schon das Rauschen zwischen den Peaks.
Diese Peaks muss ich doch jetzt irgendwie nutzen können oder?
Uns wurde gesagt dies seien die Koeffizienten der Sinusschwingungen in 
der jeweiligen Frequenz(Abszisse).
Wie kann mir die ifft Funktion genau dabei helfen?

[Sven B. (scummos)]

fourier schrieb:
> Mein Ziel ist es das Rauschen aus dem Signal zu bekommen um das
> periodische Signal zu bekommen.
> Man sieht an meiner Abbildung ja schon das Rauschen zwischen den Peaks.
> Diese Peaks muss ich doch jetzt irgendwie nutzen können oder?
> Uns wurde gesagt dies seien die Koeffizienten der Sinusschwingungen in
> der jeweiligen Frequenz(Abszisse).
> Wie kann mir die ifft Funktion genau dabei helfen?

Was willst du denn mit dem periodischen (Zeit)signal machen? Das sieht 
doch mit Sicherheit schon super aus, auch ohne Filter.

"Ich will das Signal ohne Rauschen" ist kein Ziel. Ein Ziel ist "ich 
will die Zeitpunkte der Nulldurchgänge" oder "ich will die Amplitude" 
oder sowas.

[fourier (Gast)]

Im Anhang liegt das Originalsignal.


Mein Ziel ist es mithilfe der fft die Harmonischen zu erkennen.
Das ist mir ja auch schon gelungen, nur frage ich mich ob das auch 
wirklich die Amplituden der jeweiligen Harmonischen sind. Und wie groß 
der Wert für die Frequenz der jeweiligen Harmonischen ist, welche der 
Koeffizient f
für z.B. A*sin(2*pi*f)  ist.
A ist die Amplitude aus dem Spektrum.
Im Prinzip geht es um die Synthetisierung durch Summierung der 
Sinus-Schwingungen oder man nimmt den Weg über die Funktion ifft.
Aber das bekomme ich irgenwie nicht hin.

Ich habe zum Beispiel gelesen, dass man die Amplituden der Frequenzen, 
welche keine Peaks sind Null setzt und die Peaks so lässt wie sie sind 
und dann eine
ifft macht und davon den realteil plottet .
So kann man das synthetisierte Signal darstellen. So wurde uns das 
gezeigt.

[Sven B. (scummos)]

Also irgendwelche sinnvollen Zeitsignale zu produzieren durch Nullen von 
FFT-Bins kannst du m.E. aufgeben. Du machst immer die Phasenbeziehung 
total kaputt, und zerstörst das (für die Rekonstrukstion wichtige) 
Bleeding neben den Bins, etc.

Das funktioniert genau dann so wie man es sich vorstellt, wenn alle 
enthaltenen Frequenzen genau in das gesampelte Zeitfenster reinpassen. 
Dann ist alles schön. Andernfalls ist es komplizierter.

Entweder du benutzt einen Zeitbereichs-Filter wie Mittelwert, oder 
IIR/FIR-Filtersynthese-Tools.

[fourier (Gast)]

Ich darf keine Synthesetools benutzen.

Dann bleibt mir nur die Synthetisierung durch Summierung der Sinus 
Schwingungen.
Nur was sind denn die Amplituden der einzelnen Sinus Schwingungen.
Sind das die Amplituden aus dem Betragsspektrum?

[Burkhard K. (buks)]

Sven B. schrieb:
> Also irgendwelche sinnvollen Zeitsignale zu produzieren durch Nullen von
> FFT-Bins kannst du m.E. aufgeben.

Man muss ja nicht unbedingt Nullen, wenn Absenken hilft: 
https://wiki.audacityteam.org/wiki/How_Audacity_Noise_Reduction_Works.

[Markus B. (russenbaer)]

fourier schrieb:
> Nur was sind denn die Amplituden der einzelnen Sinus Schwingungen.
> Sind das die Amplituden aus dem Betragsspektrum?

"Ja" - deshalb unter Anführungszeichen weil es auf Deine FFT 
Implementierung ankommt - im Detail wo der Faktor 1/N appliziert wird - 
bei der FFT oder bei der IFFT oder als 1/sqrt(N) bei beiden (wobei N die 
FFT Größe ist).

Ein pragmatischer Weg herauszufinden welchen Faktor Du benötigst 
(zwischen Betragsspektrum und Amplitude) ist das Du Dir einen Sinus 
erzeugst (gleiche Länge wie Dein Signal das Du analysieren willst) - mit 
bekannter Frequenz und Amplitude - und Deine FFT anwendest. Dann kannst 
Du dir direkt anschauen wie dieser Sinus im Frequenzbereich abgebildet 
wurde und dir den Faktor berechnen. Verwende für die Sinusfrequenz eine 
bin-Frequenz.
Der andere Weg führt über Dein Lehrbuch und die MATLAB Hilfe um zu sehen 
wie MATLAB die FFT implementiert hat und welche Faktoren angewendet 
werden.

Nach Deinem FFT-Bild oben kannst Du die sehr niedrigen bins schon 
Nullen. Das  ist nur Rauschen. Da machst Du Dir keine Phasenbeziehungen 
kaputt.
Was Du nicht machen darfst, und was Poster scummos mit Bleeding 
bezeichnet, und wo er Recht hat ist, das Du Dir Phasenbeziehungen 
zerstörst wenn Du von einer "starken" bin-Gruppe nur einen  bin übrig 
lässt.
Wenn Du in Deinen FFT-Plot reinzoomst wirst Du sehen das dort in so 
einer Gruppe nicht genau nur ein bin sehr "stark" ist sondern auch die 
bins links und rechts daneben (und eventuell deren Nachbarn usw. 
[abhängig von der Frequenz und von Deiner Fensterfunktion]). Die werden 
auch benötigt für die IFFT.
Du fensterst Deine Daten mit einem Rechteck-Fenster im Zeitbereich (d.h. 
Du wendest kein Fenster bei der FFT/IFFT an = Rechteckfenster). Das ist 
gleichbedeutend mit einer Faltung im Frequenzbereich mit sin(x)/x. 
Dadurch verschmierst Du etwas Energie im Bildbereich/Spektrum.


lg
Markus

[Forist (Gast)]

fourier schrieb:
> Im Anhang liegt ein JPG mit dem Amplitudenspektrum bis 10 KHz.

Warum ein JPG. Für Liniengraphik ist PNG prädestiniert.
Der JPEG-Kompressionsalgorithmus vermatscht jede scharfe Kante.