Filter Messdaten

WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

chunu

Forum-Anfänger


	Beiträge: 10

	Anmeldedatum: 20.05.15

	Wohnort: ---

	Version: ---

Verfasst am: 20.05.2015, 16:02 Titel: Filter Messdaten

Hallo,

Ich nutze den folgenden Code, um das Geschwindigkeitssignal einer Strömungsmessung zu filtern.

Code:

N = length (signal);
dF = Fs/N;
%Frequenzanalyse
signal_fft = fftshift(fft(signal))/N;
f = (-Fs/2:dF:Fs/2-dF);
% Bandpassfilter
BPF = ((f_HP < abs(f)) & (abs(f) < f_TP))';
signal_fft = BPF.*signal_fft;
% Rücktrafo
signal_gefiltert=ifft(ifftshift(signal_fft))*N;

Funktion ohne Link?

Wenn ich als signal alle Messdaten übergebe, funktioniert die Filterung. Sobald ich jedoch nur ein Teil der Daten (z.B. Sample 1 bis 2000) übergebe, bekomme ich komplexe Zahlen mit signifikant hohem Imaginärteil als Ergebnis der Filterung. Im Spektrum äußert sich dies z.B. durch ein starkes Ansteigen im eigentlich gefilterten Teil.

Habt Ihr eine Idee woran dies liegen kann, bzw. Tipps wie ich die Filterung besser umsetzen sollte?

Beste Grüße

DSP

Forum-Meister


	Beiträge: 2.117

	Anmeldedatum: 28.02.11

	Wohnort: ---

	Version: R2014b

Verfasst am: 20.05.2015, 17:01 Titel:

Ein komplexes Signal nach der Rücktransformation kann dadurch verursacht werden, wenn der Vektor signal_fft nicht exakt konjugiert symmetrisch ist.

Möglicherweise kommt es bei der FFT und der anschließenden Filterung (Faltung) zu Rundungsfehlern, wodurch die Symmetrie nicht mehr vorliegt. Das ist aber nur eine Vermutung. Evtl. spielt auch Leakage eine Rolle, da ja durch die Segmentierung möglicherweise keine kompletten Perioden im Messfenster liegen.

Eine Möglichkeit wäre folgenden Parameter zu nutzen:

Code:

y = ifft(..., 'symmetric')

Funktion ohne Link?

Im Übrigen könntest du dir den Faktor 1/N auch sparen, so fern der Code nur zur Faltung genutzt wird und nicht noch das Frequenzspektrum dargstellt wird.

Allerdings führst die Faltung auch falsch durch, nämlich eine zyklische statt einer linearen Faltung. Daher wundert es mich, dass du überhaupt in einem Fall ein richtiges Ergebnis erhältst. Beim Faltungssatz eines Signals der Länge m mit einer Filterimpulsantwort h der Länge N (nichts anderes stellt eine Filterung ja dar), ergibt sich für den Ergebnisvektor folgende Länge:

Code:

signal_filt[m+n-1] = signal[m] * h[n] % Faltungsatz

Funktion ohne Link?

Um die Faltung nun im Frequenzbereich durchführen zu können, muss unbedingt das signal wie auch h VOR der Transformation auf die Länge m+n-1 durch Anhängen von Nullen erweitert werden. Nur so ergibt sich eine lineare Faltung, die im Zeitbereich das exakt gleich Ergebnis liefert wie

Code:

signal_filt = conv(signal,h) % Faltung im Zeitbereich

Funktion ohne Link?

Siehe folgender Code zur Faltung im Freq.bereich:

Code:

function output = FFT_Faltung(sig1, sig2)
% Die Funktion führt eine schnelle Faltung mittels FFT aus
% Input:
% sig1 = Messsignal
% sig2 = Filter-Impulsantwort
% Output:
% output = gefiltertes Messignal
%--------------------------------------------------------------------------

% kopieren der Eingangsvektoren
sig1 = double(sig1(:));
sig2 = double(sig2(:));

% Faltungssatz sig(m) * h(n) = output(m+n-1):
outlength = length(sig1)+length(sig2)-1;
% nächste Zweierpotenz für FFT
fftsize = 2^nextpow2(outlength);

% Messsignal
%mit Nullen auf fftsize auffüllen
sig1 = [sig1; zeros(fftsize-length(sig1),1)];
% Berechnung des Frequenzspektrums
sig1 = fft(sig1,fftsize);

% Filter - Impulsantwort
% mit Nullen auf fftsize auffüllen
sig2 = [sig2; zeros(fftsize-length(sig2),1)];
% Berechnung des Frequenzspektrums
sig2 = fft(sig2,fftsize);

% Faltung durch Multiplikation der Spektren
% Rücktransformation in den Zeitbereich
conv_raw = ifft(sig1.*sig2);

% Rückgabe Ergebnis
output = transpose(conv_raw(1:outlength));

Funktion ohne Link?

chunu

Themenstarter

Forum-Anfänger


	Beiträge: 10

	Anmeldedatum: 20.05.15

	Wohnort: ---

	Version: ---

Verfasst am: 21.05.2015, 09:42 Titel:

Vielen Dank für Deine Antwort.
Ich werde mich insgesamt zu diesem Thema wohl noch etwas einlesen müssen. Hast Du hierzu evtl. eine Literaturempfehlung?

Mit dem 'symmetric' Parameter bei ifft scheinen die Probleme nicht aufzutreten. Jedoch hast Du recht, dass die Signale nach der Filterung nicht mehr übereinstimmen. Die Spektren von gefiltertem und ungefiltertem Signal überscheiden sich mehr als nur durch die Filterung.

Würdet Ihr mir eher zur Filterung mittels 'filter'-Funktion raten?

DSP

Forum-Meister


	Beiträge: 2.117

	Anmeldedatum: 28.02.11

	Wohnort: ---

	Version: R2014b

Verfasst am: 21.05.2015, 09:56 Titel:

chunu hat Folgendes geschrieben:

Jedoch hast Du recht, dass die Signale nach der Filterung nicht mehr übereinstimmen. Die Spektren von gefiltertem und ungefiltertem Signal überscheiden sich mehr als nur durch die Filterung.

Was meinst du damit genau? Ein Bild wäre zur Erklärung sicherlich hilfreich.

Zwar in Englisch aber auf jeden Fall für Anfänger sehr geeignet: http://www.dspguide.com/ch18/2.htm

Die Kapitel 9 und 18 beschäftigen sich mit der Faltung. Ansonsten z.B. das Buch Digitale Signalverarbetiung von Kammeyer/Kroschel. Oder Google benutzen...Stichtworte FFT Faltung, Overlap-Add Methode

Ob du nun filter/conv oder die Faltung im Frequenzbereich durchführst. Alle sollten das gleiche Ergebnis liefern. Für die filter(B,A,signal) benötiges du halt die Filterkoeffizenten A und B der Polynome. Oder du erstellt dir eine Impulsantwort des Filters. Dann kannst du entweder conv nutzen und über FFT. Je nach Länge des Signals und Filters hat die FFT Faltung eine deutlich geringe Ausführungszeit. Daher wählt man ja auch den Umweg über die Transformation, da sich eine einfache Multiplikation ergibt.

chunu

Themenstarter

Forum-Anfänger


	Beiträge: 10

	Anmeldedatum: 20.05.15

	Wohnort: ---

	Version: ---

Verfasst am: 21.05.2015, 10:18 Titel:

Ein Bild ist angehängt. Die blaue Kurve wurde mit folgendem Code erzeugt. Rot ist das ursprüngliche Signal

Code:

N = length (y);
dF = Fs/N;

y_fft = fftshift(fft(y));
f = (-Fs/2:dF:Fs/2-dF);

BPF = ((f_HP < abs(f)) & (abs(f) < f_TP))';
y_fft_filt = BPF.*y_fft;

y_filt=ifft(ifftshift(y_fft_filt),'symmetric');

Funktion ohne Link?

VGL Filter.png

Beschreibung:

Download

Dateiname:

VGL Filter.png

Dateigröße:

18.92 KB

Heruntergeladen:

332 mal

DSP

Forum-Meister


	Beiträge: 2.117

	Anmeldedatum: 28.02.11

	Wohnort: ---

	Version: R2014b

Verfasst am: 21.05.2015, 10:29 Titel:

Aber wie passt das Bild zu deiner Aussage?

Zitat:

Die Spektren von gefiltertem und ungefiltertem Signal überscheiden sich mehr als nur durch die Filterung.

Da ich nicht weiß welchen Bereich genau der Bandpass durchlassen soll (f_HP, f_TP ??) kann ich nur raten. Für micht sieht es eher so aus als würde nur der Hochpass wirken, bzw. die Dämpfung des Tiefpasses nur marginal ist. Ansonsten erkenne ich da jetzt erstmal nichts Ungewöhnliches.

Wie kommst du eigentlich darauf das die Y-Achse der Leistungsdichte entspricht?

chunu

Themenstarter

Forum-Anfänger


	Beiträge: 10

	Anmeldedatum: 20.05.15

	Wohnort: ---

	Version: ---

Verfasst am: 21.05.2015, 13:03 Titel:

Ich meinte mit der Aussage, dass sich die spektren mehr als nur durch den gefilterten Bereich unterscheiden.

Sry. Hab natürlich weitere Infos vergessen. f_HP = 1, f_TP = samplingrate. Soll primär nur als Hochpass dienen, um z.B. das Signal vom langsammen Drift zu befreien.

Habe das Spektrum schnell mit pwelch erzeugt und deswegen Leistungsdichte genannt.

DSP

Forum-Meister


	Beiträge: 2.117

	Anmeldedatum: 28.02.11

	Wohnort: ---

	Version: R2014b

Verfasst am: 21.05.2015, 13:12 Titel:

Ok, jetzt habe ich deine Aussage verstanden.

Möglich ist aber auch, dass die Unterschiede durch die falsche Faltung verursacht werden. Eine korrekten Faltungs-Algorithmus habe ich dir ja bereits gepostet. Dort musst du ja nur sig2 durch dein Filter ersetzen. Oder in deinem Code die Erweiterung von Nullen für y hinzufügen und entsprechend auch N anpassen damit der Frequenzvektor wieder passt.

Zitat:

Habe das Spektrum schnell mit pwelch erzeugt und deswegen Leistungsdichte genannt.

Ok, dann ist das Label richtig Wink

chunu

Themenstarter

Forum-Anfänger


	Beiträge: 10

	Anmeldedatum: 20.05.15

	Wohnort: ---

	Version: ---

Verfasst am: 21.05.2015, 14:02 Titel:

Habe Deinen Code wie beschrieben angewendet. Ergebnis im Anhang.
Irgendetwas stimmt nicht Smile

VGL Filter2.png

Beschreibung:

Download

Dateiname:

VGL Filter2.png

Dateigröße:

16.52 KB

Heruntergeladen:

305 mal

Signal_Output.png

Beschreibung:

Download

Dateiname:

Signal_Output.png

Dateigröße:

9.68 KB

Heruntergeladen:

281 mal

DSP

Forum-Meister


	Beiträge: 2.117

	Anmeldedatum: 28.02.11

	Wohnort: ---

	Version: R2014b

Verfasst am: 21.05.2015, 14:13 Titel:

Bitte auch mal den Code dazu posten.

chunu

Themenstarter

Forum-Anfänger


	Beiträge: 10

	Anmeldedatum: 20.05.15

	Wohnort: ---

	Version: ---

Verfasst am: 21.05.2015, 14:33 Titel:

Vielen Dank für Deine schnellen Antworten und die Hilfe!

Code:

y = u(nfirst:nlast); %Segment aus Gesamtsignal

N = length (y);
dF = Fs/N;

f = (-Fs/2:dF:Fs/2-dF);

BPF = ((f_HP < abs(f)) & (abs(f) < f_TP))';

sig1 = y;
sig2 = double(BPF(:));

% Faltungssatz sig(m) * h(n) = output(m+n-1):
outlength = length(sig1)+length(sig2)-1;
% nächste Zweierpotenz für FFT
fftsize = 2^nextpow2(outlength);

% Messsignal
%mit Nullen auf fftsize auffüllen
sig1 = [sig1; zeros(fftsize-length(sig1),1)];
% Berechnung des Frequenzspektrums
sig1 = fft(sig1,fftsize);

% Filter - Impulsantwort
% mit Nullen auf fftsize auffüllen
sig2 = [sig2; zeros(fftsize-length(sig2),1)];
% Berechnung des Frequenzspektrums
sig2 = fft(sig2,fftsize);

% Faltung durch Multiplikation der Spektren
% Rücktransformation in den Zeitbereich
conv_raw = ifft(sig1.*sig2);

% Rückgabe Ergebnis
output = transpose(conv_raw(1:outlength));

Funktion ohne Link?

DSP

Forum-Meister


	Beiträge: 2.117

	Anmeldedatum: 28.02.11

	Wohnort: ---

	Version: R2014b

Verfasst am: 21.05.2015, 15:04 Titel:

Ist ungetestet, also k.A. ob es funktioniert. Du machst aber schon den Fehler das dein Frequenzvektor nicht zu sig1 passt, da dieser anders angeordnet ist. Daher war es bei dir schon richtig ffftshift zu nutzen. Außerdem ist dein Filter doch schon im Frequenzbereich, bzw. es gewichtet eben das Messignal im Freq.bereich. Wozu also BPF noch transformieren?

Code:

y = u(nfirst:nlast); %Segment aus Gesamtsignal

sig1 = y;

% Faltungssatz sig(m) * h(n) = output(m+n-1):
outlength = length(sig1)+length(sig2)-1;
% nächste Zweierpotenz für FFT
fftsize = 2^nextpow2(outlength);

% Messsignal
% mit Nullen auf fftsize auffüllen
sig1 = [sig1; zeros(fftsize-length(sig1),1)];
% Berechnung des Frequenzspektrums
sig1 = fftshift(fft(sig1,fftsize));

% Filter - Frequenzbereich
dF = Fs/fftsize;

f = (-Fs/2:dF:Fs/2-dF);

BPF = ((f_HP < abs(f)) & (abs(f) < f_TP))';
sig2 = double(BPF(:));

% Faltung durch Multiplikation der Spektren
% Rücktransformation in den Zeitbereich
conv_raw = ifft(ifftshift(sig1.*sig2));

% Rückgabe Ergebnis
output = transpose(conv_raw(1:outlength));

Funktion ohne Link?

Mittlerweile würde ich aber eher über das Filter an sich nochmals nachdenken. Einen Bandpass kannst du auch mittels Matlabfunktionen und Tools erstellen. Es gibt z.B. das "fdatool". Dann kannst du einfach die filter() Funktion verwenden.

http://de.mathworks.com/help/dsp/re.....?searchHighlight=bandpass

chunu

Themenstarter

Forum-Anfänger


	Beiträge: 10

	Anmeldedatum: 20.05.15

	Wohnort: ---

	Version: ---

Verfasst am: 21.05.2015, 16:11 Titel:

Hier eine Umsetzung mit 'butter' und das zugehörige Ergebnis. Die Dämpfung reicht mir jedoch nicht aus. Sobald ich jedoch die Ordung über 7 erhöhe, kommt nichts Sinnvolles (Spektrum liegt um einige Größenordungen höher) heraus.
Im Bereich um 10 Hz stimmt doch auch etwas nicht oder?

Vieleicht eher grundsätzlich die Frage was Ihr mir vorschlagen würdet, um mein Signal vom Trend zu befreien.

Code:

y = u(nfirst:nlast)
f_hp_n = 1/(0.5*Fs);

[b a] = butter(5,f_hp_n,'high');
output = filter(b,a,y);

Funktion ohne Link?

Filter.png

Beschreibung:

Download

Dateiname:

Filter.png

Dateigröße:

18.13 KB

Heruntergeladen:

261 mal

DSP

Forum-Meister


	Beiträge: 2.117

	Anmeldedatum: 28.02.11

	Wohnort: ---

	Version: R2014b

Verfasst am: 21.05.2015, 16:26 Titel:

Da ich gerade kein Matlab zur Verfügung habe, kann ich es leider selber nicht testen. Vermutlich ist die geringe Dämpfung auf die sehr geringe Grenzfrequenz von 1Hz zurückzuführen. Mit höherer Ordnung wird das Filter dann instabil, weshalb kein sinnvolles Signal mehr herauskommt.

Evtl. mal die Grenzfrequenz etwas anheben. Mal mit 10Hz beginnen und sich die Dämpfung anschauen. Einfacher ist dazu aber das fdatool, da man dort gleich den Frequenzgang des Filters darstellen kann und sieht welche Dämpfung es bei den eingestellten Parametern erreicht.

Ansonsten mal statt butter ein Chebyshev Filter testen:

Code:

[b,a] = cheby2(n,Rs,Ws,ftype)

Funktion ohne Link?

Oder weg von den IIR Filtern zu den FIR. Die können nicht instabil werden. Da gäbe es z.B. die sogenannten Window-Sinc Filter. Ebenfalls im fdatool realisierbar oder hier mal schauen:

http://www.gomatlab.de/window-sinc-filter-t19156.html

chunu

Themenstarter

Forum-Anfänger


	Beiträge: 10

	Anmeldedatum: 20.05.15

	Wohnort: ---

	Version: ---

Verfasst am: 22.05.2015, 16:05 Titel:

Das Verhalten tritt bei allen getesteten Filtern auf.
Ich vermute mittlerweile dass es an der sehr geringen Grenzfrequenz von 2 Hz und der hohen Abtastrate von 2000 Hz liegt.
Jemand eine Idee, wie ich diese Parameter umsetzen kann?

	Gehe zu Seite 1, 2 Weiter
Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:	Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen

Impressum | Nutzungsbedingungen | Datenschutz | FAQ |

RSS

Hosted by:

Copyright © 2007 - 2025 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.