www.mikrocontroller.net

Forum: Digitale Signalverarbeitung / DSP autocorrelation, matlab/octave


Autor: daniel (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hallo Leute,

ich verstehe nicht warum bei meiner T=1 peridosichen
Funktion die Autokorrelationsfunktion nicht periodisch ist.

Die Andeutungen der Parabeln sind zwar zu sehen,
aber die Peaks der 100% Übereinstimmung im Abstand T
werden nicht erreicht.

Im Anhang hab ich einen Screenshot beigefügt
wie es bei mir aussieht.
Rot ist die Autokorrelation.

und hier ist mein Code
************************************
#!/usr/bin/env octave

x1 = -1:0.1:1;
x5 = repmat(x1, 1, 5);
t = linspace(0, 5, length(x5));

for i=0:length(x5)
    clf;
    hold on;
    plot(t, x5, 'b');
    plot(t, shift(x5, -i), 'g');

    plot(t, autocor(shift(x5, -i))', 'r');

    pause(0.2);
    grid;
end

pause;

************************************

Was hab ich übersehen oder Gedankenfehler gemacht?

Grüsse, Daniel

Autor: pumpkin (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Moin Daniel,

was stört dich? Die Höhe oder der Abstand delta_t?

Ich hab mal mit ein Beispiel mit Matlab gemacht:
y=sin(0:.1:6*pi);
corr=xcorr(y,y);
plot(corr)

Sieht wunderbar aus. Dass die peaks nach links und rechts kleiner werden 
ist klar, das Integral wird durch die Verschiebung "kürzer" - mal ganz 
doof ausgedrückt.

pumpkin

Autor: daniel (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hi,

ah wahrscheinlich hast du recht
wärne es nicht 5 Perioden sondern unendlich nach rechts und links
und würde man sich dann einen Ausschnitt betrachten, dann
wären die Peaks alle gleich hoch.

Ich könnte probieren mehr Perioden zu nehmen und die autocor
aber über ein paar wenige anschauen

Gruss, Daniel

Autor: daniel (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
.. und im Code

octave:6> x1=-1:0.1:1;
octave:7> x1000=repmat(x1, 1, 1000);
octave:8> xc=autocor(x1000)';
octave:9> t = linspace(0,1,length(x1)*5);
octave:10> plot(t,x1000(1:length(t)));
octave:11> hold on;
octave:12> grid
octave:13> plot(t,xc(1:length(t)),'r');

Autor: pumpkin (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
So ist es. Das ganze sieht dann so aus:
y=sin(0:0.1:60*pi);
xcorr(y,y);
plot(ans(1697:2072))

pumpkin

Autor: pumpkin (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
:^)   Zwei Doofe, ein Gedanke.

Autor: daniel (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
wenigstens ist jetzt klarer ;)

Also im Grunde handelt es sich bei Matlab und Octave mit
den Funktionen autocor/xcor um Autokorrelation der Energiesignale.
Ausserhalb des gegebenen x Intervalles wird Funktion zu 0 angenommen
und Integral von -inf bis +inf berechnet.
Bin gerade beim Lernen und versuch die ganzen unterschiedlichen
Definitionen voneinander abzugrenzen.
Sowas wie Autokorrelation für Leistungssignale gibt's
wahrscheinlich nicht in Matlab?

Autor: pumpkin (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Richtig! Bisher habe ich keine für Leistungssignale entdeckt (ich hab 
aber auch noch nicht danach gesucht). Kann man sich aber selbst 
schreiben, es gibt ne nette Funktion namens circshift() die ich da 
verwenden würde. Kannst du dir ja mal zu Gemüte ziehen.

pumpkin

Autor: Andreas Schwarz (andreas) (Admin) Benutzerseite Flattr this
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Die periodische Autokorrelation entspricht der inversen DFT vom 
Leistungsdichtespektrum, also ifft(abs(fft(x1)).^2).

EDIT: hatte abs() vergessen

Autor: pumpkin (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Stimmt, Wiener-Theorem? Schon wieder verdrängt...

pumpkin

Autor: Matthias (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
xcorr(x1,x1) liefert übrigens auch die zyklische Autokorrelation. Für 
eine "normale" Autokorrelation muss man noch genügend (ich glaub 
length(x1)-1) 0en an jedes Signal anhängen.

Autor: Andreas Schwarz (andreas) (Admin) Benutzerseite Flattr this
Datum:

Bewertung
0 lesenswert
nicht lesenswert
xcorr macht eine normale Autokorrelation:
> x=[3 2 1];
> xcorr(x)
ans =

   3   8  14   8   3

Zyklische Autokorrelation:
> ifft(abs(fft(x)) .^ 2)
ans =

  14  11  11

Herleiten kann man sich das mit:
Faltung(x, y) = IDFT(DFT(x) * DFT(y))
Korrelation(x, y) = Faltung(x, reverse(y)) = IDFT(DFT(x) * 
DFT(reverse(y)))
Autokorrelation(x) = Faltung(x, reverse(x)) = IDFT(DFT(x) * 
DFT(reverse(x))) = IDFT(DFT(x) * conj(DFT(x))) = IDFT(abs(DFT(x)).^2)

Die autocor-Funktion arbeitet mit autocov() und ist mit etwas Vorsicht 
zu genießen:
> autocor(x)
ans =

   1.00000
   0.00000
  -0.50000

EDIT: ein reverse() vergessen

Autor: daniel (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
hab gerade nachgedacht, es gilt nur für x reel oder?

kleiner Fehler Andreas?
Autokorrelation(x) = Faltung(x, -x)

Aber Danke für diesen Trick!
Gruss, Daniel

Autor: Andreas Schwarz (andreas) (Admin) Benutzerseite Flattr this
Datum:

Bewertung
0 lesenswert
nicht lesenswert
y muss eigentlich noch konjugiert werden, dann funktioniert das auch mit 
komplexen Zahlen:

Korrelation(x, y) = Faltung(x, conj(reverse(y)))

Den Fehler habe ich oben korrigiert.

Autor: Shadravan Fontanov (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hallo,

hier ein Thread mit dem aehnlichen Thema:
http://www.mathworks.com/matlabcentral/newsreader/...

Gruesse

 Shadravan

Autor: Peter K. (peter26)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hi ...

Also bezüglich der Einhüllenden -> Rückgang der Amplitude zu den Rändern 
hin gibts 2 Möglichkeiten dies zu reduzieren:

.) biased cross correlation -> hier steigt jedoch die Varianz des 
Schätzers der Korrelation nach aussen an. (am besten einfach ausprobiern 
-> Matlab-Normalisierung)

.) frequeny based correlation -> FFT^-1(FFT(w) . *FFT(w)) fällt der 
Abfall nach aussen weg da ja eine periodische Fortsetzung angenommen 
wird.

lg Peter

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.