FIR Filter um Sinus zu detektieren

Gast #2643577
Lesenswert?

In Matlab gibt es eine Funktion mit dem Namen goertzel. Leider besitze 
ich diese Funktion nicht. Nun habe ich mir gedacht ich programmiere den 
Goertzlalgorithmus mit der filter Funktion nach.

Siehe Übertragungsfunktion:
http://www.mathworks.de/help/toolbox/signal/ref/goertzel.html
1
%N = 256
2
for i=0:1:N-1
3
  f = filter([1 -exp(-i*2*pi*i/N)],[1 -2*cos(2*pi*i/N) 1],1:N);
4
end;

Müsste diese Implementierung so stimmen ?
Wie füge ich nun das Sinussignal (100kHz mit N = 256 Abtastwerte) hinzu 
?
Gast #2644112
Lesenswert?

Sehe ich nicht, wie das gehen soll, was macht die Filterfunktion?

Musst Du nicht einfach nur für jeden Wert die jeweils beiden 
zurückliegenden Werte aufgreifen und mit dem COS-Faktor multiplizieren, 
wie es in WIki steht? Ich meine das englische Wiki!
Gast #2644174
Lesenswert?

Tip: Etwas generischer geht's mit der Signal Processing toolbox so:
1
fgoer = tf([0.5], [ 1, -(r*kcos1), r*r ]);
2
filtered = filter(tf.num, tf.den, data);

tf plus entsprechende Argumente baut die Transfer-Funktion für den 
Goertzel.
r ist typischerweise fast 1 (0.9999...), kcos1 der cos(2*pi*k/N)-Term 
für die gewünschte Frequenz.
Gast #2644645
Lesenswert?

Hallo Strubi,

danke für deine Hilfe. Ich nicht nachvollziehen wie du auf diese Zeile 
kommst:
>>fgoer = tf([0.5], [ 1, -(r*kcos1), r*r ]);

Mir fehlt eine genaue Herleitung des Goerzel Algorithmus.
Gast #2644804
Lesenswert?

Ich habe den Filter so implementiert wie auf der Internetseite 
http://www.google.de/url?sa=t&rct=j&q=sa2011-lab3-fourier&source=web&cd=1&ved=0CCcQFjAA&url=http%3A%2F%2Fwww.emt.tugraz.at%2Fsystem%2Ffiles%2FSA2011-Lab3-Fourier.pdf&ei=8y6RT8aELo2GhQfkhb2fBA&usg=AFQjCNEok5d4o0Efmwge58wWGKVdjsM_fw 
beschrieben (Seite 9).
1
f=filter([1],[1 -exp(-j*2*pi*k/N)],y)

y sind meine Eingangsdaten
N ist Konstant
und k ist der Laufindex von k=0... N-1

Nach der Berechnung mache ich dann eine Betragsbildung mit dem Befehl 
abs.
Stimmt meine Implemntierung in Matlab ?
Gast #2644889
Lesenswert?

Techniker schrieb:
> Stimmt meine Implemntierung in Matlab ?

Wie wäre es wenn Du Dir einfach ein Testsignal erzeugst, z.B. so was:
1
f=50;            %Signalfrequenz
2
fs=100*f;        %Abtastfrequenz
3
t=0:1/fs:1/f;    %Zeitvektor erstellen
4

5
y= sin(2*pi*f*t)+sin(2*pi*f*t)+sin(4*pi*f*t);

Dann lässt Du den Filter drüber laufen.
Anschließend ne FFT und Du weist es.

Grüße
Gast #2644972
Lesenswert?

Irgendwie funktioniert das ganze nicht so richtig.
Was mache ich da noch falsch ?
1
f=50;            %Signalfrequenz
2
fs=100*f;        %Abtastfrequenz
3
t=0:1/fs:1/f;    %Zeitvektor erstellen
4

5
y= sin(2*pi*f*t);
6
max_y = max(abs(y))*1.1; 
7
subplot(2,1,1);
8
fig = figure(1); 
9
plot(t,y);
10
axis([0 t(end) -max_y max_y]);
11
grid on;
12

13
N = 128;      % Abtastwerte
14
k = 0:1:N-1;
15
g=filter([1],[1 -exp(-j*2*pi*k/N)],y);
16

17
H = fft(g,N); 
18
amplH = abs(H);
19
amplitudengang = fftshift(amplH/N);
20
subplot(2,1,2);
21
stem(amplitudengang);
22
grid on;
Persönliche Seite #2645592
Lesenswert?

Mathegenie schrieb:
> Was funktioniert nicht richtig? Siehst du Deine Frequenz nicht, oder
> siehst Du gar nichts?

Kann es sein, dass Du Deine Frequenz deshalb nicht findest, weil Du zu 
grob abtastest?

Die Schrittweite der Testfrequenzen mit denen Du rechnest, muss in einem 
sinnvollen Verhältnis zur Samplezahl stehen, weil mit zunehmender 
Samplezahl die DFT-Rechnung immer selektiver wird und Dir so eine 
Frequenz durchrutscht. Ich habe im Anhang ein altes Beispiel hinkopiert:

Die Schrittweite von 8 führt dazu, dass bei der blauen Kurve mit den 
Punkten die Frequenz einigermassen getroffen wird und sich ein peak 
bildet, bei 10 Schritten liegt man haarscharf daneben und der peak geht 
im Rauschen unter. Bei der gröber auflösenden, türkisen Kurve wurde die 
Berechung schon nach 2048 von 4096 samples abgebrochen.

Umgekehrt führt die Überlegung dahin, dass Du Deine Sanplezahl nicht zu 
gross wählen darfst, wenn Du Deine Frequenz nicht genau triffst, oder Du 
bekommst eben das Ergebnis, daß das Matching gering ist.


Techniker schrieb:
> Würde es eigentlich Sinn machen
> einen Bandpass FIR-Filter zu verwenden ?
Der Bandpass ist nicht anderes als eine Summe von Frequenzen als Maske, 
die gleichzeitig auf den Datenstrom losgelassen werden. Mathematisch ist 
der Bandpass ein Hochpass, von dem ein TP abgzogen wurde. In den 
Koeffizienten zeigt sich das durch ein Integral der oberen Grenzfrequenz 
hin zu 0, von der der Koeefizientensatz unterhalb der GF weggelassen 
wurde.

Wenn Du dir einander annäherst, hast Du am Ende eine Frequenz wie oben.

Wenn Du Deine Frequenz nicht genau kennst, ist der BP gfs besser.
Angehängte Dateien:
Gast #2648342
Lesenswert?

Hallo Jürgen,

danke für deinen Beitrag. Du meinst ich taste zu grob ab. Auf einem 
System kann habe ich nun mal pro Sinusperiode 128 Abtastwerte. Das 
müsste doch reichen. Ich habe die Vermutung das irgendwo noch ein Fehler 
drin ist. Bin leider noch nicht weitergekommen.
Gast #2648344
Lesenswert?

Mit dem Goertzelalgorithmus möchte ich erkennen, ob der Sinus zum 
Beispiel mit 100Khz vorhanden ist oder nicht. Wenn Impulsstörungen oder 
auch Rauschen additiv mit dem Sinus überlagert sind, soll trotzdem der 
Sinus von 100Khz detektiert werden.
Persönliche Seite #2701128
Lesenswert?

Techniker schrieb:
> danke für deinen Beitrag. Du meinst ich taste zu grob ab. Auf einem
> System kann habe ich nun mal pro Sinusperiode 128 Abtastwerte.

Ich denke, wenn dann diesbezüglich "zu fein". Je mehr Punkte je Sinus, 
desto selektiver wird diese Maske, mit der Du den Datenstrom ansiehst. 
Wenn diese Punktezahl steigt, muss auch die "Maskenanzahl" steigen.

In der Grafik habe ich das veranschaulicht:

* Blau: Sehr genau abgetastet, man erkennt den Verlauf des Signals
* Rot: Mit einem 7-tel der Aufösung abtastet: Das Maximum rutscht durch
* Grün: Mit einem 7-tel der Auflösung bei 7-facher Breite
* Hellgrün genau abgetastet mit 7facher Breite

Der grüne ist also grob/breit genug, um das Signal noch zu erfassen und 
fein genug, um den Verlauf noch erahnen zu können.

Blau wäre etwa die obere Grenze. Grün etwa die untere Grenze.
.
Angehängte Dateien:

Antwort schreiben

Bitte melde dich an, um einen Beitrag zu schreiben.

oder

Mit Google-Account einloggen

Die Registrierung ist kostenlos und dauert nur eine Minute.

Jetzt registrieren