Verfasst am: 31.05.2012, 16:04
Titel: FFT und IFFT wo ist der Hund begraben ?
Hallo Community,
Ich tüftel grade mit der FFT und der IFFT rum und nun passiert folgendes :
Code Version 1
Code:
%Signal Fourier Transformieren
Fs = 1; %Sampling Frequenz in Hz
T = 1/Fs; %Sample Time
L = 151; %Länge des Signals
t = (0:L-1)*T; %Zeit Vector
NFFT = 2^nextpow2(L); %Next Pow of 2 from Length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
%Signal Fourier Transformieren
Fs = 1; %Sampling Frequenz in Hz
T = 1/Fs; %Sample Time
L = 151; %Länge des Signals
t = (0:L-1)*T; %Zeit Vector
NFFT = 2^nextpow2(L); %Next Pow of 2 from Length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
Du hast 151 Werte und übergibst der FFT aber mit NFFT die nächsthöhere 2er Potenz = 256....fft(y,NFFT). Somit hängt die Funktion automatisch vor der Transformation an y entsprechend die fehlenden Werte als NULLEN an, um auf die Länge NFFT zu kommen. Bei der Rücktransformation werden die Nullen dann naturlich wieder hergestellt. Diese müsstest du entweder abschneiden, oder NFFT nicht übergeben. Ist L dann keine 2er Potenz wird der langsame DFT Algorithmus verwendet.
Edit: Schau dir mal in der Skriptecke das "umfassende FFT Bsp." an...deine Skalierung stimmt nämlich nicht ganz (Gleichsignalanteil und bei Fs/2)
Erstmal vielen Dank für die Hilfe
Auf die sache mit dem NFFT hätte ich auch kommen können steht ja auch im manual.
Auch das vorgeschlagene Script hilft schon ein wenig. Das Problem ist nun, das ich das fft Signal gerne mit nem Butterwort filter Filtern will aber das müsst ich ja auch das oben stehende "f" anwenden sonst geht das ja nicht. Aber ich kann das "f" ja nicht mit der ifft zurück transformieren
Wo ist denn in dem oberen Code ein Filter? f ist doch nur deine Frequenzachse.
Es gibt mehrere Möglichkeiten zum Filtern...
- Funktion filter() benutzen
- Faltung im Zeitbereich von Signal und Filter-Impulsantwort mittels Funktion conv()
-Faltung im Freq.-bereich...Filter-Impulsantwort mittels FFT transformieren und dann mit Y[f] multiplizieren.
Zu jeder Möglichkeit gibt es hier im Forum genügend Bsp...einfach mal suchen
nochmals danke für die Antworten !
Ich formuliere die sache mit dem Filter etwas genauer:
Code:
%Signal Fourier Transformieren
Fs = 1; %Sampling Frequenz in Hz
T = 1/Fs; %Sample Time
L = 151; %Länge des Signals
t = (0:L-1)*T; %Zeit Vector
NFFT = 2^nextpow2(L); %Next Pow of 2 from Length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
Do = 20; %Grenzfrequenz
n = 3; %Grad des Filter
P = 1; %nachlesen wegen 2 Dimensionen und so
for ind=1:1:151
grad3(ind,1)=butterworthfilt(Do,n,ind,P);
end figure(2);
plot(grad3,'green');
hold on;
n = 2;
for ind=1:1:151
grad2(ind,1)=butterworthfilt(Do,n,ind,P);
end figure(2);
plot(grad2,'red');
n = 1;
for ind=1:1:151
grad1(ind,1)=butterworthfilt(Do,n,ind,P);
end figure(2);
plot(grad1,'blue');
Ohne Daten kann ich deinen Code nicht testen. In grad ist dann auch wirklich die Impulsantwort des Filters gespeichert? Diese brauchst du nämlich, um die Faltung durchzuführen. Du machst auch noch einen Fehler bei der Faltung im Freq.bereich. Bei der Faltung von x[n] und h[m] kommt am Ende eine Länge von n+m-1 heraus. Beachtet man dies nicht, wird eine zyklische statt der linearen Faltung durchgeführt und damit ist das Ergebnis falsch. Du musst x[n] und h[m] vor der Transformation auf die Länge n+m-1 mit Nullen erweitern. Siehe folgende Funktion für die Faltung...
Code:
function output = FFT_Faltung(sig1, sig2) % Die Funktion führt eine schnelle Faltung mittels FFT aus % Input: % sig1 = Messsignal % sig2 = Filterimpulsantwort % Output: % output = gefiltertes Messignal
%--------------------------------------------------------------------------
% kopieren der Eingangsvektoren
sig1 = double(sig1(:));
sig2 = double(sig2(:));
% 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);
% Ausgabe
output = transpose(conv_raw(1:outlength));
Führe doch erstmal die Faltung richtig durch und schaue dir das Ergebnis an. Für die Zukunft wäre ich/wir auch dankbar, wenn du an deinen Plots eine Achsenbeschriftung hättest. So muss man immer raten, was hier eigentlich dargestellt ist.
Ja ich probiere das mal aus und schreibe dann nochmal.
Entschuldige die Sache wegen den Plots das ist etwas lieblos.
Aber du hast schon recht, ich werde das ab jetzt berücksichtigen. Ich kann ja nicht erwarten das man sich die Mühe macht mir zu antworten, aber mir selber nicht die Mühe machen eine ordentliche Frage zu stellen.
Ich habe dein Programm jetzt mal getestet...was soll denn grad sein? Die Filter-Impulsantwort ist es jedenfalls nicht. Es sieht in etwa wie eine Frequenzantwort eines Filters aus. Aber warum solltest du sie dann noch mal mittels FFT transformieren. Ehrlich gesagt kann ich auch mit deiner Formel in der Butterworth Fkt. nichts anfangen.
So ich habe das Jetzt noch einmal etwas "verschönert"
Ein Paar Hintergrund Informationen Vorweg
Mein gemessenes "Signal" sind Ortungspunkte in einem Raum also zum Zeitpunkt t war ich an einer einer coordinate (x,y). Meiner Meinung nach lässt sich dies nicht als 2 Dimensinal betrachten, denn so würde ich ja aus einer weißen fläche mein eigentliches Signal raus filtern denn x,y ist nur eine Coordinate ohne einen weiteren Wert oder Ähnliches.
Butterworth
der filter sieht etwas komisch aus da ich ihn aus einem buch für Bildverarbeitung habe und eine Dimension weggestrichen habe und ein paar Überreste geblieben sind. Gängigerweise gibt es ja dieses vorgehen mit dem verdoppeln der Bildgröße und dem vertauschen der Ecken, damit das Spektrum zentriert wird.
Natürlich kann ich den Vorgefertigten Matlab filter des SignalProc. Toolkits nehmen, aber da ich dies unter anderem für mein Studium mache finde ich das das schlichte anwenden der Funktion keinen mehrwert fürs verständnis mit sich bringt.
[edit]Die Filterfunktion ist augenscheinlich ja auch richtig, zumindest optisch ist es gleich dem was in der Literatur zu finden ist [/edit]
Code:
%Datei Laden
%Datei Pfad
xlsread('C:\Users\Grossmann_HK2\Downloads\Rundweg1.xls');
x = ans(:,11); %x Coordinaten aus Datei holen
y = ans(:,12); %y Coordinaten aus Datei holen
%Signal Fourier Transformieren
Fs = 1; %Sampling Frequenz in Hz
T = 1/Fs; %Sample Time
L = 151; %Länge des Signals
t = (0:L-1)*T; %Zeit Vector
NFFT = 2^nextpow2(L); %Next Pow of 2 from Length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
Do = 20; %Grenzfrequenz
n = 3; %Grad des Filter
P = 1; %nachlesen wegen 2 Dimensionen und so
for ind=1:1:151
grad3(ind,1)=1/(1+((ind-(P/2))/Do)^(2*n));
end subplot(3,3,4);
plot(grad3,'green');
xlabel('');
ylabel('');
title('');
hold on;
n = 2;
for ind=1:1:151
grad2(ind,1)=1/(1+((ind-(P/2))/Do)^(2*n));
end
Mit Bildverarbeitung kenne ich mich ehrlich gesagt nicht so aus.
Da es ja aber hier ein 1D Signal ist, sind meine vorherigen Aussagen über die Faltung im Freq.-bereich richtig. Leider beherzigst du diese Ratschläge nicht, weshalb ich nicht großartig Lust habe hier noch weiterzuhelfen. Der Faltungssatz sind Grundlagen der Mathematik und in dieser Form sicherlich auch der Signalverarbeitung. Nachlesbar in vielen Lektüren zu diesem Thema.
Ansonsten sehe ich hier keinen Fehler bei der Skalierung, da ja weder bei dem Signal = Y[f] noch bei der Filterfreq.-antwort grad1,2,3 = H[f] eine Skalierung stattfindet. Diese kann man sich bei der Faltung ja auch sparen. Somit kann die Verstärkung, die in den gefilterten Signalen zu sehen ist, nur vom Filter herkommen.
Ok...ich habe dir ja noch einen Ansatz genannt, wo in meinen Augen das Problem liegt. Ob zyklische oder die richtige lineare Faltung hat mit deinem Skalierungsproblem nämlich nichts zu tun.
Wie kommst du eigentlich darauf, dass deine Abtastfreq. 1 Hz beträgt? Gleichzeitig willst du das Signal dann bei 20 Hz Tiefpass filtern. Das kann ja nicht gehen.
Einstellungen und Berechtigungen
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