DFT mit Matlab berechnen

[Roman (Gast)]

Hallo,
ich verzweifel hier schon seit ein paar Stunden und zwar habe ich 
Messdaten in einer .txt file. In diesem File sind zwei Spalten mit 
jeweils einen Zeitpunkt und wiederrum ein dazugehöriger Amplitudenwert. 
Diese möchte ich nun mit Hilfe der DFT Transformieren.
Leider klappt es vorne und hinten nicht. Ich bekomme nichts sinnvolles 
aus meinem Graphen raus.
Das Eingangssignal ist ein Sinussignal bei f = 400kHz und der Amplitude 
5V

clear all; close all; clc; 
A = importdata('sinus_test2.txt', ';'); % die Daten sind mit einem Semikolon getrennt 
B = importdata('sinus_test2.txt', ';'); 


tmin =10.e-6;           % Mit tmin und tmax soll später eine Fensterung durchgeführt werden 
tmax =100e-6;             
tmin_counter = 0; 
tmax_counter = 0; 

 for i=1:length(A)  % ermittelt den Bereich, in dem t min und tmax liegt 
    if(A(i,1)<tmin) 
      tmin_counter = tmin_counter +1; 
    end 
    if(A(i,1)>tmax) 
      tmax_counter = tmax_counter +1; 
    end 
 end 
tmax_counter = length(A)-tmax_counter; 
  
 for i=1:tmax_counter-tmin_counter% entnehme die werte aus Datei zwischen tmin und tmax 
    t(i)=B(i+tmin_counter,1); 
    Pin(i)=B(i+tmin_counter,2); 
    x(i)=A(i+tmin_counter,2); 
 end 
%---------------------------------- 
fs = 1/(t(2)-t(1)); % das entspricht meinen Abtestwerten 

%x=sin(2*pi*100e3*t); ALternativ habe ich mal hier probiert einen sinus zu plotten 
plot(t,x); 
hold on 
% Sample frequency (Hz) 
m = length(x);          % Window length 
n = pow2(nextpow2(m));  % Transform length 
y = fft(x,n);           % DFT 
f = (0:n-1)*(fs/n);     % Frequency range 
power = y.*conj(y)/n;   % Power of the DFT 
plot(f,power) 
xlabel('Frequency (Hz)') 
ylabel('Power') 
title('{\bf Periodogram}') 
plot(f(1:floor(n/2)),power(1:floor(n/2)))

Ich denke mal das ich einen Fehler bei der Frequenznormierung habe... 
aber ich weiß nicht wie ich das anders llösen kann,

Und ja ich kenne die folgenden Seiten schon
http://www.gomatlab.de/fft-umfassen.....l-t777,highlight,dft.html
http://de.mathworks.com/help/matlab.....ourier-transform-fft.html

[Lukas K. (carrotindustries)]

Dass deine Zeitwerte nicht äquidistant sind, ist dir schon aufgefallen?

[Guido C. (guidoanalog)]

Hallo,

entspricht die erste Spalte in "sinus_test.txt" der Zeit in Sekunden?
Falls ja ergibt sich aus den ersten zwei Abtastzeitpunkten eine
Abtastfrequenz von:
f = 1/(1.00014687248443E-5 s - 1.00006066506084E-5 s)
f = 1.1600e+09 Hz = 1,16 GHz
Tastest Du wirklich so schnell ab?

Des Weiteren sind Deine Abtastwerte nicht äquidistant (vergl. Anhang).

Last but nut least: Vermeide in Matlab Schleifen wann immer es geht, die
Programmlaufzeit wird es Dir danken.

Mit freundlichen Grüßen
Guido

[Roman (Gast)]

Lukas K. schrieb:
> Dass deine Zeitwerte nicht äquidistant sind, ist dir schon
> aufgefallen?

Ehrlich gesagt, nicht das ist echt peinlich und ich habe mich schon die 
ganze Zeit gefragt wieso ich beim plotten der Funktion so einen 
komischen murks raus bekomme...

Scheiße, jetzt muss ich mir überlegen wie ich das problem umgehen kann. 
WÜrdet ihr empfehlen, dass ich die Funktion mit Matlab interpolieren 
soll und dann mit der FFt transformiere?

Guido C. schrieb:
> f = 1/(1.00014687248443E-5 s - 1.00006066506084E-5 s)
> f = 1.1600e+09 Hz = 1,16 GHz
> Tastest Du wirklich so schnell ab?

Ja das stimmt soweit. Die Messergebnisse kommen von einer 
Simulationssoftware die leider keine FFT/DFT Funktionen besitzt. 
Deswegen kommt man auch auf den extrem hohen Abtastwert von über 1 GHz.

[Lukas K. (carrotindustries)]

Roman schrieb:
> Scheiße, jetzt muss ich mir überlegen wie ich das problem umgehen kann.
> WÜrdet ihr empfehlen, dass ich die Funktion mit Matlab interpolieren
> soll und dann mit der FFt transformiere?

Mir war da gerade mal fad...
https://gist.github.com/carrotIndustries/095ef01df03cd1e59512
Ist zwar python, die Ideen aber dieselben.

[Guido C. (guidoanalog)]

Hallo Roman,

anbei mein Denkanstoß.

'[MATLAB]
clear all; close all; clc; 
A = importdata('sinus_test.txt', ';'); % die Daten sind mit einem Semikolon getrennt 

Ax=linspace(A(1,1),A(end,1),length(A(:,1)))';
Ay=interp1(A(:,1),A(:,2),Ax);

figure;
plot(Ax,Ay);
grid on;
hold on;

figure;
f_A = 1/(Ax(2)-Ax(1));
Ax_FFT = f_A/length(A(:,1))*(0:length(A(:,1))-1)';
Ay_FFT = abs(fft(Ay));
semilogy(Ax_FFT,Ay_FFT);
grid on;
hold on;
'[/MATLAB]

Mit freundlichen Grüßen
Guido

[Roman (Gast)]

Recht herzlichen Dank für die Hilfe. Das hat mir schon einiges 
geholfen... Das einzige was ich noch nicht genau hinbekommen habe ist 
die Amplitudeninformation. Wie kann ich diese genau Dimensionieren? 
Normalerweise würde man durch die Länge der Mesppunkte -1  Dividieren... 
aber ich bekomme dann bei meinem SInussignal nichts sinnvolles raus....

[Lukas K. (carrotindustries)]

Roman schrieb:
> Das einzige was ich noch nicht genau hinbekommen habe ist
> die Amplitudeninformation.

Bei einer reellen FFT musst du durch n/2 teilen. Wenn du davor noch 
fensterst, musst du die Abschwächung des Fensters (Mittelwert des 
Fensters) mit berücksichtigen. Zudem kann die Wahl des Fensters die 
Amplitudengenauigkeit erheblich beeinflussen. Wenn du genaue Amplituden 
willst, nimm ein Flattop-Fenster. Stell dir vor, dein Peak fällt genau 
zwischen zwei FFT-bins. Wenn das Fenster im Frequenzbereich bereits bei 
halber bin-breite abgefallen ist, hast du einen Amplitudenfehler. Das 
Flattop-Fenster ist im Bereich um ein halben Bin - wie der Name sagt - 
flach.

[Roman (Gast)]

Wenn ich den oben genannten Code verwende, reicht es da nicht aus, wenn 
ich die X-Achse skaliere und dann durch einen Faktor dividiere ?

[Guido C. (guidoanalog)]

Hallo,

Roman schrieb:
> Wenn ich den oben genannten Code verwende, reicht es da nicht aus, wenn
> ich die X-Achse skaliere und dann durch einen Faktor dividiere ?

auf welchen Code beziehst Du Dich? Den Code von Lukas oder den Code von 
mir? Falls Du Dich auf meinen Code beziehst, die x-Achse sollte dort 
bereits richtig skaliert sein.

Mit freundlichen Grüßen
Guido

[Roman (Gast)]

Guido C. schrieb:
> '[MATLAB]
> clear all; close all; clc;
> A = importdata('sinus_test.txt', ';'); % die Daten sind mit einem
> Semikolon getrennt
>
> Ax=linspace(A(1,1),A(end,1),length(A(:,1)))';
> Ay=interp1(A(:,1),A(:,2),Ax);
>
> figure;
> plot(Ax,Ay);
> grid on;
> hold on;
>
> figure;
> f_A = 1/(Ax(2)-Ax(1));
> Ax_FFT = f_A/length(A(:,1))*(0:length(A(:,1))-1)';
> Ay_FFT = abs(fft(Ay));
> semilogy(Ax_FFT,Ay_FFT);
> grid on;
> hold on;
> '[/MATLAB]

Ich beziehe mich auf deinen, Guidos, Code.
Die X-Achse ist soweit auch richtig skaliert. Nun möchte ich gerne 
wissen wuie ich das mit der Fensterung machen soll ? Warum reicht es bei 
deinem Code nicht aus, einfach den den Ay Wert durch einen bestimmten 
Faktor zu teilen ?

Z.b. Als Refernsignal kann ich ja z.B. eine SInusamplitude nehmen mit 
Peak Wert von 1 W -> daraus Mittelwert Bestimmen -> Danach weiß man ja 
was rauskommen sollte...

[Schneider007 (Gast)]

Hallo Roman,

Windowfunctions nimmt man doch, damit man bessere (eindeutigere) Peaks 
bekommt. Da brauch ich das Signal auch nur einfach mit der 
Windowfunction multiplizieren, dass die der Anfang und das Ende des 
Signals nicht so abrupt starten bzw. enden. Ich würde übrigens noch die 
erste Zeile streichen.

clear all; close all; clc; erschwert das debuggen unnötig und alte 
Fehlermeldungen werden gelöscht.

[Schneider007 (Gast)]

deshalb würde ich die erste Zeile streichen.
Ansonsten viel Erfolg.
mfg

[Robin (Gast)]

Ok, habe ich getan, aber leider passt die Leistung immernoch nicht

Bx=linspace(B(1,1),B(end,1),length(B(:,1)))';
By = 15 *sin(2*pi*1e6*Bx); % 10W -> 38,49dbm , 15W-> -> 40,255 dBm
f_B = 1/(Bx(2)-Bx(1));      % Abtastzeitpunkte festlegen

w = flattopwin(length(Bx));
wsignal = By(:).*w; % Fensterung vornehmen
ws=sum(w);
Y= 10*log10(1000*abs(fft(wsignal)/ws)); % logarithmieren
Bx_FFT = f_B/length(B(:,1))*(0:length(B(:,1))-1)';
plot(Bx_FFT,Y);

Nun bekomme ich ein Peak bei meiner 1 MHz raus aber mit einer 
Ausgangsleistung von 38,75 dBm. Da fehlen ja rein theoretisch immer noch 
fast 2 dB.

[Schneider007 (Gast)]

Matlab hat auch schon vorgefertigte Windowfunctions, ne. Hamming oder so 
zB, musste mal googeln. Wenn du hier keine Hilfe findest, würde ich dir 
gomatlab.de empfehlen. Da war ich auch ne Zeit unterwegs, da sind auch 
Profis die helfen. Ich glaube aber, wenn hier kein Fehler bei der 
Wichtungsfunktion ist, dann ist ein Rechenfehler und dazu kenne ich mit 
fft (die du ja hier machst und nicht dft) zu wenig aus. Ob sich einer 
auf gomatlab.de damit genug auskennt weiß ich auch nicht, aber versuchen 
kann man es ja. Viel Erfolg.

[Roman (Gast)]

Ja, vielen Dank schon einmal!