Verfasst am: 23.07.2009, 16:58
Titel: schnell und einfach einen Sinus fitten?! Frequenz gesucht!
Hallo,
ich suche nach einer Möglichkeit die mit automatisch die Frequenz von Interferenzstreifen ausgibt.
Mit den fit Möglichkeiten von Matlab muss ich aber immer zu viele input Parameter vorgeben und sind diese nicht richtig gewählt wirds garnix (siehe Bild, grün sind die Datenpunkte!) Es gibt ja genügend Software (Origin und Co) welche sowas schnell und brauchbar fitten kann bzw. mir die Frequenz ausgeben kann, sodass ich dies in mein Programm integrieren kann ohne die Frequenz da immer von Hand grob einzugeben.
Ts = diff(zeit(1:2)); % Abtastzeit
N = length(zeit); % Länge des Datenvektors
f = [0:floor((N-1)/2)] / (N*Ts); % Frequenzvektor für Spektrum-Plot % Fouriertransformierte
SIGNAL = fft(signal); % Fouriertransformation
SIGNAL = SIGNAL / N; % Normierung
SIGNAL = [SIGNAL(1)2*SIGNAL(2:floor((N-1) / 2) + 1)]; % begrenzen auf f_max
ergebnis = [SIGNAL; f];
Ziel ist eine diskrete FFT eines vorliegenden, zeitlichen Signalverlaufes. Die Funktion benötigt einen Signal und einen Zeitvektor und erstellt daraus die Fouriertransformierte des Signals und einen Frequenzvektor. Mit den Befehlszeilen
Sich das Frequenzspektrum anzuschauen ist zwar eine gute Methode, wie ich davon aber jetzt die Frequenz geliefert bekomme verstehe ich noch nicht. Gibt es einen peakfinder bzw. am besten nen Gauss-fit der automatisch die dominante Frequenz (ist ja nur ein peak) aus dem Spektrum jetzt bestimmen könnte?
Das ganze soll und muss ja automatisiert ablaufen, da ich sehr viele Datensätze mit immer völlig verschiedener Frequenz habe!
Du hast das Problem aber auch schon benannt. Mit dieser Methode erhalte ich eine Periode von 79 samplen für das Beispiel oben, wobei das gut 10% Fehler hat! Manuell gemessen kommt man da auf 70-72 sampel pro Periode.
Da es sich bei den Datensätzen um Interferenzstreifen auf einem Bild handelt kann ich die Sampleanzahl natürlich nicht erhöhen. Wenn ich das richtig sehe, hat sich die Fouriertranformation für das Problem damit erledigt. Die Frequenzen der Interferenzstreifen sind immer sehr gering, von daher kann ich so wohl keinen genauen Wert erhalten.
Vielleicht hab ich das Problem falsch umschrieben: Man nehme ein Bild von Streifen. Wie ist der Abstand der Streifen (in pixeln)?
Natürlich komme ich auch über die Frequenz drauf, aber wie gesagt, FFT schachtelt die Frequenz viel zu grob. Ein Fehler von 10% ist zuviel, da man es ja gar manuell viel genauer hinkriegt.
ich hab mal eben auch mit der fit funktion rumgespielt und muß leider bestätigen, dass man die Frequenz bei so einem Sinus schon sehr gut kennen muss um ein ordentliches fitting zu bekommen.
Vielleicht kann man ja mit dem FFT-Verfahren die Frequenz als Startwert für ein fitting bekommen
als Startwert für die Amplitude würde ich das maximum des Signals vorgeben
und der Startwert der Phasenlage ist eh recht egal... also 0
das Problem ist, man kennt die Frequenz und Amplitude garnicht! Man weiß nur, dass es ein Sinus^2 ist, aber alles andere ändert sich von Datensatz zu Datensatz! Den Verschiebungswinkel kennt man auch nicht, der variiert wie alles andere!
Aufgrund der Amplitudenschwankungen bringt ein Standard-Fit nicht viel und ein Polynomfit hat natürlich nicht mehr eine einzelne Frequenz, so ganz einfach wirds wohl nicht zu lösen sein. Die beste Idee, die ich bisher hab ist einfach die Extremstellen-Abstände zu mitteln, schön ist das nicht
figure('Position',[10,100,1200,800]) hold on
plot(k,signal,'r') plot(k,signal_neu,'b') plot(Nu.c,Nu.value,'mO','LineWidth',2,'MarkerSize',10) hold off
grid on
Problem dabei ist, dass durch Messrauschen gerne mal aus einem Nulldurchgang eine Reihe von Nulldurchgängen wird... die müßte man aber doch erkennen können indem man ein Fenster vorgibt, in dem nur EIN Nulldurchgang sein kann...
kommen wir der Lösung mit dem Ansatz evtl näher?
anzahl = length(Nu.value)
fenster = 5; % Anzahl Abtastwerte innerhalb der nur EINE Nullstelle liegen darf
delta = diff(Nu.c); % Abstände der Nullstellen ermitteln
count = 0;
i = 1;
Nu_neu.c = Nu.c;
while i<=length(delta)% zu häufige Nullstellen rauskicken if delta(i) > fenster
delta(i) = delta(i) + round(count/2); % versuch die verschmierte Nullstelle zu mitteln
Nu_neu.c(i) = Nu_neu.c(i) - round(count/2);
count = 0;
end if delta(i) <= fenster
delta(i) = [];
Nu_neu.c(i) = [];
i = i - 1;
count = count + 1;
end
i = i + 1;
end
figure('Position',[10,100,1200,800]) hold on
plot(k,signal,'r') plot(k,signal_neu,'b') plot(Nu.c,Nu.value,'mO','LineWidth',2,'MarkerSize',10) plot(Nu_neu.c,Nu_neu.value(),'kX','LineWidth',2,'MarkerSize',10) hold off
grid on
ich hoffe man kann das nachvollziehen... ist jetz nur mal so eine Schnell-Lösung... aber mit den vorgegebenen Daten klappts ganz OK...
die Variable 'delta' gibt jetzt die Abstände der Nulldurchgänge aus...
je nach Fenster-Breite dürfte das ganz ordentlich funktionieren...
Jetzt musst Du nur noch die Abtastrate wissen und kannst die gesuchte Frequenz berechnen.
Die Auswertung der Maxima hat den Nachteil, dass Maxima von Sinus-Wellen recht flach aussehen... da wird die Bestimmung immer etwas unscharf... Die Nullstellen-Bestimmung ist da besser, WENN man ein mittelwertfreies Signal hätte!
die Mittelwertberechnung mit mean ist nur richtig, wenn du auch komplette Sinuswellen in den Messdaten hast! Wenn du natürlich (wie im vorliegenden Fall) die Sinusschwingung abschneidest... hmm dann wirds wieder etwas falsch... da hilft dann nur MEHR aufnehmen, so dass der Fehler klein wird. Also ein längeres Meßintervall...
irgendwie so würde ich versuchen der Lösung näher zu kommen
EINE letzte ganz andere Idee fällt mir da noch ein...
wird nicht die Autokorrelation eingesetzt um in verrauschtesten Signalen irgendeine Periodizität zu finden...
evtl kann man durch diese Herangehensweise punkten ?
habe ich mir bisher allerdings noch nicht überlegt...
die Idee mit den Nulldurchgängen ist gut, hatte ich auch, funktioniert aber nicht, da das Signal nicht Symmetrisch sein muss!
Ich hab am Wochenende folgendes gebastelt was ganz gut funktioniert. Ich schau mir einfach die Abstände der Extrema an. Ne gute statistische Analyse mach ich da aber noch nicht, hab aus der Hilfe nicht direkt ersehen können wie ich z.B. ne Standardabweichung auf einen Mittelwert ausgeben kann....
Code:
%% Evaluate fringe frequency % ============================== clc;clear all;close all;
% input DATA
b=imread('images/luft.bmp');
sortierbreite=10; %amount of lines on top and bottom to consider
image_scale=1; %compared to original image(important for scale/phase only)
G=double(b(:,:,2)); %Farbe Grün in Zeile x Spalte x 2
G=G/(max(max(G))); %normalizes to 1;
maxima_difference=[0];
minima_difference=[0];
%% Drop peripheral extrema % first and last extremum are likely not real ones due to cutoff!
if minima(1,1)<maxima(1,1);
minima=minima(:,[2:size(minima,2)]);
else
maxima=maxima(:,[2:size(maxima,2)]);
end;
if minima(1,size(minima,2))>maxima(1,size(maxima,2));
minima=minima(:,[1:size(minima,2)-1]);
else
maxima=maxima(:,[1:size(maxima,2)-1]);
end;
%% Calculate difference
ifsize(minima,2)>1; %if just one extrememum left, no distance can be calculated
for i=1:size(minima,2)-1;
minima_difference(h,i)=minima(1,i+1)-minima(1,i);
end;
end;
ifsize(maxima,2)>1;
for i=1:size(maxima,2)-1;
maxima_difference(h,i)=maxima(1,i+1)-maxima(1,i);
end;
end;
%=============================================
%=============================================
% SAME FOR LOWER EDGE
%=============================================
%=============================================
%for h=size(G,1)-sortierbreite+1:size(G,1);
for h=1:sortierbreite;
%% Drop peripheral extrema % first and last extremum are likely not real ones due to cutoff!
if minima(1,1)<maxima(1,1);
minima=minima(:,[2:size(minima,2)]);
else
maxima=maxima(:,[2:size(maxima,2)]);
end;
if minima(1,size(minima,2))>maxima(1,size(maxima,2));
minima=minima(:,[1:size(minima,2)-1]);
else
maxima=maxima(:,[1:size(maxima,2)-1]);
end;
%% Calculate distance of extrema
ifsize(minima,2)>1; %if just one extrememum left, no distance can be calculated
for i=1:size(minima,2)-1;
minima_difference(sortierbreite+h,i)=minima(1,i+1)-minima(1,i);
end;
end;
ifsize(maxima,2)>1;
for i=1:size(maxima,2)-1;
maxima_difference(sortierbreite+h,i)=maxima(1,i+1)-maxima(1,i);
end;
end;
%% Plotting Origianl & Smoothened Data + Extrema subplot(3,1,2),plot(signal,'black.'),title('Lower edge'); hold on;
subplot(3,1,2),plot(signal2,'blue');
subplot(3,1,2),plot(minima(1,:),minima(2,:),'r*') subplot(3,1,2),plot(maxima(1,:),maxima(2,:),'g*')
%average maxima without zeros!!
counter=0;
maxima_difference_mean=0;
for s=1:size(maxima_difference,1);
for u=1:size(maxima_difference,2);
if maxima_difference(s,u)~=0;
maxima_difference_mean=maxima_difference_mean+maxima_difference(s,u);
counter=counter+1;
elseif maxima_difference(s,u)==0;
maxima_difference(s,u)=NaN;
end;
end;
end;
maxima_difference_mean=maxima_difference_mean/counter;
%average minima without zeros!!
counter=0;
minima_difference_mean=0;
for s=1:size(minima_difference,1);
for u=1:size(minima_difference,2);
if minima_difference(s,u)~=0;
minima_difference_mean=minima_difference_mean+minima_difference(s,u);
counter=counter+1;
elseif minima_difference(s,u)==0;
minima_difference(s,u)=NaN;
end;
end;
end;
minima_difference_mean=minima_difference_mean/counter;
nein, das ist keine gute Idee, da die Amplitude moduliert sein kann und es oft ist!
habe vergessen zu erwähnen, dass ich die function "peakdet" von Eli Billauer verwende! (umbenannt in "tools_peakdetection".
http://www.billauer.co.il/peakdet.html
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