Verfasst am: 27.09.2012, 15:01
Titel: Filterung und FFT passen nicht immer zusammen
hallo an alle wissende
ich bin hier schon ewig am rumprobieren und habe einen neuen rückschlag, bei dem ich alleine nicht weiter komme...
also ich hab ein zeitsignal, welches ich durch eine FFT mit fenster transformiere und anschließend noch ein BP anwenden will. das funktionierete auch bisher. umindest solange ich es nur auf einen kleinen signalausschnitt angewendet habe( L= 6000) bei deutlich längeren signalen (L=250000) haut dann mein bandpass nciht mehr hin, weil die Vektorlängen nicht emhr überein stimmen.
sicher hat es was damit zu tun, dass man ja für die fft mit
N=2^nextpow(L)
arbeitet, so dass sich die vektorlänge des amplitudengangs immer ändert. aber ich hab jetzt keinen ansatz, wie ich das auf meine zeit-signale umwandle. irgendwie muss sich dass ja für den bandpass auch anpassen lassen... will das programm später in einer gui anwenden, sowohl, wenn ich nur einen ausschnitt betrachte, als auch für das gesamte zeitsignal...
hier der code
Code:
% Laden des Files
%-----------------------------------
filename = uigetfile('*.xls','Auswahl der gewünschten Zeitreihe');
[nums, txt] =xlsread(filename);
%----------------------------------
%Zeitbereich
% ----------------------------------
x =nums(:,1); % x-Vektor
y =nums(:,2); % y-vektor
Ta = 6.65e-6; % Abtastrate
fa = 1/Ta; % Abtastfrequenz
fn = fa/2; % Nyquistfrequenz
L = length(x); % Länge des Zeitsignals
N = 2^nextpow2(L); % gewünschte FFT-Länge (N=2^x, sonst wird der DFT-Algorithmus verwendet!)
df = fa/N; % Frequenzauflösung
%---Darstellung Datensatz ---
% max. Amplitude zur Skalierung der graphischen Darstellung feststellen:
max_y = max(abs(y))*1.1;
fig = figure(1); %fig 1 plot(x,y);
%axis([0 N -max_y max_y]) title('Zeitreihe') ylabel('Spannung') xlabel('Zeit') grid
% Anhang an die bereits erfolgte Untersuchung % -------------------------------------------
win = hann(L)';
y_win = y'.*win; % Fensterung ohne Amplitudenkorrektur
%y_win = y'.*win*N/sum(win); % Fensterung mit Amplitudenkorrektur
max_y = max(abs(y_win))*1.1;
% Berechnung der FFT % ------------------
H = fft(y_win, N);
% Berechnung des Amplitudengangs aus dem komplexen Frequenzvektor H:
amplH = abs(H);
% Darstellung des interessierenden Frequenzbereichs des Amplitudengangs (0...fn) und % daran angepasste Amplitudenskalierung (Normierung auf N/2):
%-------------------------------------------------------------------------------
amplitudengang = [amplH(1)/N amplH(2:N/2)/(N/2)]; % DC-Bin auf N normieren!
title('Amplitudengang nach Fensterung') ylabel('Amplitude') xlabel(['Auflösung: ',num2str(df),' Hz Frequenz in Hz']) grid
%-----------------------------
% Bandpass
%-----------------------------
%Ordnung n:
n = 10000; % Anzahl der Nullstellen in der kompl. Ebene
pass1 = 1000; % -6 dB Grenzfrequenz in Hz
pass2 = 28000; % -6 dB Grenzfrequenz in Hz
Wn = [pass1/fn pass2/fn];
FIRkoeff = fir1(n, Wn, 'bandpass'); % BP-Filter
% Zunächst den Nennerkoeffizientenvektor erstellen: % (nötig zur korrekten Rechnung mit freqz, grpdelay, zplane)
a = zeros(size(FIRkoeff)); % oder:
%a = zeros(1, n+1);
a(1) = 1;
% freqz über Aufruf mit Koeffizientenvektoren
cntW = 4096;
% Standardmäßig (ohne zusätzliche Parameterangabe cntW) werden nur 512 Frequenzpunkte % berechnet! [Hfir, Wfir] = freqz(FIRkoeff, a, cntW); % mit 0 <= Wfir <= pi
%df = Wfir(2)/pi*fn; % oder: df = fn/length(Wfir); % df (delta_f) ist der Abstand zwischen den Frequenzpunkten in Hz % d.h. df entspricht der Frequenzauflösung in Hz
amplFIR = abs(Hfir);
phaseFIR = angle(Hfir);
% Darstellung des Amplitudengangs (linear):
fig = figure(fig+1); %fig
plot(x_fn, amplFIR, 'b');
%plot(Wfir/pi*fn, amplFIR, 'b');
axis([0 fn -0.11.1]);
title('Amplitudengang des FIR-Filters') ylabel('Verstärkung') xlabel(['Auflösung: ',num2str(df),' Hz Frequenz in Hz']) hold on;
plot(Wn*fn,1,'rx') hold off;
%multiplikation mit frequenzgang
y_BP = amplitudengang'.*amplFIR;
Das Problem liegt vermutlich in der Faltung des Signals mit der Impulsantwort des Filters. Es ist richtig, dass diese im Freq.-bereich in eine Multiplikation übergeht. Allerdings muss man dabei eine zyklische Faltung vermeiden, da sonst Fehler entstehen. Dazu müssen die Längen von Signal und Impulsantw. gleich sein. Schau mal hier bei der Funktion Faltung...
Hab überlegt, das Problem mit zero padding zu umgehen. Also meinen Bandpass einfach mit Nullen am Ende aufzufüleln, damit die Signale gleich lang sind.
Wenn tmp1 z.B. 1000 ist, erzeugt "zeros(tmp1)" eine [1000x1000] Matrix. Das wird mit etwa 8MB Speicherplatz noch nichgt den Speicher sprengen, aber ich vermute, Dein tmp1 ist deutlich größer, und Du wolltest eigentlich einen Vektor erstellen, also:
Mal abgesehen von diesem Fehler ist das zeropadding so immer noch nicht korrekt, bzw. die Faltung liefert dir wieder ein falsches Ergebnis.
Du musst den Faltungssatz beachten...
Aus dem Signal x[n] der Länge n und der Imulsantwort h[m] der Länge m wird ein Ausgangssignal y der Länge n + m - 1!
Daher müssen vor der Transformation die beiden Signal x und h auf die Länge n+m-1 durch zeropadding gebracht werden. Siehe die Funktion FFT_Faltung.m aus dem Link. Nur so wird eine zyklische Faltung vermieden. Der Nachteil dieser Methode ist, dass dadurch eben viele Nullen angefügt werden müssen, vor allem wenn n>>m ist. Hierzu gibt es dann gesonderte Verfahren (Overlap Add/Save) wodurch das Zeropadding reduziert werden kann, da das Signal x in Segmente geteilt wird. Aber hier sind dann Korrekturen an den Enden der Segmente notwendig um das richtige Ergebnis zu erhalten.
Pepsi
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 01.10.2012, 16:08
Titel:
ah, ok danke! ich werd mich die tage mal ransetzen und schauen, dass ich das hinbekomme! vielen dank schonmal
Pepsi
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 02.10.2012, 15:02
Titel:
also die variante mit dem zero-padding hat sich in diesem fall wohl erledigt...
durch die deutlich größere länge des signals verringert sich die frequenzauflösung. da sich amplFIR jedoch nicht ändert (was ich nicht ganz verstehe, weil ich sich fn ja anpassen müsste...?), kommt es zur verschiebung der grenzfrequenzen im amplitudengang meines bp-filters (zumindest, so wie es oben programmiert ist)...
da kommt man wohl um die faltung nciht drum rum...
Das kann nicht sein...da bist du auf dem Holzweg. Du darfst das fn eben nicht ändern...das ändert sich nämlich auch gar nicht. Die Abtastfrequenz fs bleibt ja immer noch gleich.
Die Freq.-auflösung df = fs/N (N= Anzahl Messwerte)
nur diese ändert sich durch Zeropadding.
Nutze die Funktion FFT_Faltung.m aus dem oben angegebenen Link bzw. siehe Code unten. Du kannst das Ergebnis sowohl mit den Matlab Funktionen conv wie auch filter vergleichen und wirst sehen, dass sie identsich sind (mal abgesehen von der Länge des Ergebnisvektors). Die Funktion filter schneidet den Ausschwingvorgang einfach ab...tatsächlich hat das Ergebnis aber die Länge n+m-1.
Außerdem darfst du doch nicht nur die Beträge für die Multiplikation nehmen. Du musst den Ergebnisvektor mit Real- und Imaginärteil aus der FFT für die Faltung nehmen.
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(:));
hatte die letzten tage viel anderes zu tun, aber hab mcih jetzt nochmal in ruhe rangesetzt... das mit der faltung funktioniert wunderbar, wer hätte es gedacht
für alle, die auch dran hängen hier mein code an den oben beschriebenen auszug angepasst mit den plots zur veranschaulichung:
% Laden des Files
%-----------------------------------
filename = uigetfile('*.xls','Auswahl der gewünschten Zeitreihe');
[nums, txt] =xlsread(filename);
%----------------------------------
%Zeitbereich
% ----------------------------------
x =nums(:,1); % x-Vektor %Zeitsignal
y =nums(:,2); % y-vektor
Ta = 6.65e-6; % Abtastrate
fa = 1/Ta; % Abtastfrequenz
fn = fa/2; % Nyquistfrequenz
L = length(x); % Länge des Zeitsignals
N = 2^nextpow2(L); % gewünschte FFT-Länge (N=2^x, sonst wird der DFT-Algorithmus verwendet!)
df = fa/N; % Frequenzauflösung
%---Darstellung Datensatz ---
% max. Amplitude zur Skalierung der graphischen Darstellung feststellen:
%max_y = max(abs(y))*1.1;
fig = figure(1); %fig 1 plot(x,y);
%axis([0 N -max_y max_y]) title('Zeitreihe') ylabel('Spannung') xlabel('Zeit') grid
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
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.
RSS