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

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

  
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Transfer Function mit tfestimate und manuell
 

ThKo26
Forum-Century

Forum-Century

Beiträge: 184
Anmeldedatum: 21.09.17
Wohnort: ---
Version: 2015b
     Beitrag Verfasst am: 15.06.2020, 13:24     Titel: Transfer Function mit tfestimate und manuell
  Antworten mit Zitat      
Hallo zusammen,

es gibt weider ein Thema.

Ich versuche die Ergebnisse zwischen der tfestiamte und dem manuellen Vorgang zu verstehen bzw. richtig zu machen.

mein Code:
Code:

%% TF
clc; clear; close all

analysis = 'YAWRSWA';%


        load('C:\Users\thko\Desktop\YAWRSWA.mat')    % YAW Rate
        data       = YAWR;


fs         = FS;
WindowPara = [0.1881, -0.36923, 0.28702, -0.13077, 0.02488]; % Literatur
Resolution = 0.35;                                           % Literatur
Overlap    = 0.785;                                          % Literatur

NT = length(data)* (1/ fs);         % Länge des Zeitfensters in Sekunden
Abstand = 1/NT;                     % Abstand der einzelen Linien Resolution
N = ((0.5 * fs)/ Resolution)* 2;    % Blocklänge
Res = ((0.5 * fs)*2)/ length(data); % Abstand der einzelen Linien Resolution (max. possible)
% Do some argument validation
if N > length(data)
    disp([num2str(Res), ' Hz (max. Frequenzauflösung)'])
    error(message('Frequenzauflösung passt nicht zu Blocklänge'));
end
WindowSize = N;
nfft = 2^nextpow2(length(data));
% gewünschte FFT-Länge (N=2^x, sonst wird der DFT-Algorithmus verwendet!)
% --> nicht unbeidngt ntowenidg
Overlap    = round(Overlap* WindowSize);
z          = 2*pi/WindowSize* (1:WindowSize-1); % x-Vector --> function gencoswinTK
Window     = WindowPara(1) + WindowPara(2)*cos(z)...
    + WindowPara(3)* cos(2*z) + WindowPara(4)*cos(3*z) + ...
    WindowPara(5)* cos(4*z); % Matlab Code von gencoswin Function/ BMW


       

        [magn, freq]   = tfestimate(SWA, YAWR, Window, Overlap, [], fs);
        [coh, freqcoh] = mscohere(SWA, YAWR, Window, Overlap, [], fs);
       
                maxlim = 4;
        figure;
        subplot(3,1,1); hold all; grid on
        p1 = plot(freq, real(magn));
        xlabel('Frequenz [Hz]'); ylabel('Amplitude')
        xlim([0 maxlim]); xticks(0:1:maxlim)
       
        subplot(3,1,2); hold all; grid on
        p1 = plot(freq, rad2deg(angle(magn)));
        xlabel('Frequenz [Hz]'); ylabel('Phase [deg]')
        xlim([0 maxlim]); xticks(0:1:maxlim)
       
        subplot(3,1,3); hold all; grid on
        p1 = plot(freqcoh, coh);
        xlabel('Frequenz [Hz]'); ylabel('Coherence [-]')
        xlim([0 maxlim]); xticks(0:1:maxlim)
        %         ylim([.8 1]); yticks(.8:.05:1)
       
%% ALTERNATIVE
        x = SWA;
        y = YAWR;


        window = 1;
       
        x1 = filter(ones(1,window)/window,1,x); % no filter, --> raw data
        y1 = filter(ones(1,window)/window,1,y);
%         figure;
%         subplot(1,2,1); hold all; grid on
%         plot(x); plot(x1)
%         subplot(1,2,2); hold all; grid on
%         plot(y); plot(y1)

        u =x1;
        y =y1;
       
        dt = 1/fs;
       
        N = length(u);
        U = fft(u);
        Y = fft(y);
        H1 = U .* Y ./ U.^2; % the H1 transfer function estimate
        MagData = (abs(H1));
        PhaseData = rad2deg(angle(H1));
       
        N2 = floor(N/2);
        FreqData = [0:N-1]/(N*dt); % Hz
       
        figure;
        subplot(211); plot(FreqData(1:N2),MagData(1:N2));
        hold on; grid on
        subplot(212); plot(FreqData(1:N2),PhaseData(1:N2));
        hold on; grid on
       
        % Filter der frequenzabhängigen Daten
       
        MagData   = filter(ones(1,window)/window,1,MagData);
        PhaseData = filter(ones(1,window)/window,1,PhaseData);
        FreqData  = filter(ones(1,window)/window,1,FreqData);
       
%         subplot(211); h=semilogx(FreqData(1:N2),MagData(1:N2));
%         axis([0.01 4 -inf inf])
%         subplot(212); p=semilogx(FreqData(1:N2),PhaseData(1:N2));
%         axis([0.01 4 -inf inf])
%         set(h,'Color','g');
%         set(p,'Color','g');
       
        n = length(x1);
        c = fft(x1)/n;
        amp = 2*abs(c);
        amp(1) = amp(1)/2;
        m = floor(N2/2);

        % Ermittelt die letzte Frequenz mit mind. 10% Maximalamplitude!
        frequenz=linspace(0,(m-1)*fs/n,m);
        f_end=max(frequenz);
        for i=1:m
            i=m+1-i;
            if amp(i)>0.1*max(amp)
                f_end=frequenz(i);
                break
            end
        end
       
        for i=1:N2
            if FreqData(i)>f_end
                MagData=MagData(1:i);
                PhaseData=PhaseData(1:i);
                FreqData=FreqData(1:i);
                N2=i;
                break
            end
        end
       
        low_lim=0.2; %Nicht größer als 0.2Hz setzen! Sonst fehlen natürlich die Gain_0_2Hz Werte!
        for p2=1:N2-1
            if FreqData(p2)<=low_lim && FreqData(p2+1)>low_lim
                break
            end
        end
       
        %[~,p2]=find(FreqData<0.1+0.05 & FreqData>0.1-0.05);
        %Umrechnung und abschneiden von ersten Werten (Sonst
        %extreme Skalierung!)
        MagData    = MagData(p2(1):N2);
        FreqData   = FreqData(p2(1):N2);
        PhaseData  = PhaseData(p2(1):N2);
        N2         = N2-p2(1);
        figure
        subplot(1,2,1); hold all; grid on
        h = plot(FreqData(1:N2),MagData(1:N2),'linewidth',2); hold all
        axis([low_lim maxlim -inf inf]);
        smooingdata =smoothdata(MagData, 'movmedian');
        h = plot(FreqData(1:N2),smooingdata(1:N2),'linewidth',2);
        xlabel('Frequenz [Hz]'); ylabel('Magnitude [deg]')
        subplot(1,2,2); hold all; grid on
        h=plot(FreqData(1:N2),PhaseData(1:N2),'linewidth',2); hold all
        axis([low_lim maxlim -inf inf]);
        smooingdata =smoothdata(PhaseData, 'movmedian');
        h=plot(FreqData(1:N2),smooingdata(1:N2),'linewidth',2);
        xlabel('Frequenz [Hz]'); ylabel('Phase [deg]')


 


leider sind die Ergebnisse nicht identisch.
Fragen:
wo liegt der Fehler in der Anwedung von tfestimate ?
Wieso startet die phase bei einem positiven Wert und nicht bei 0?

Bitte um Unterstützung

Viele Grüße
Thomas

YAWRSWA.mat
 Beschreibung:

Download
 Dateiname:  YAWRSWA.mat
 Dateigröße:  35.44 KB
 Heruntergeladen:  216 mal
Private Nachricht senden Benutzer-Profile anzeigen


ThKo26
Themenstarter

 Forum-Century

Forum-Century

Beiträge: 184
Anmeldedatum: 21.09.17
Wohnort: ---
Version: 2015b
     Beitrag Verfasst am: 17.06.2020, 07:15     Titel:
  Antworten mit Zitat      
Ich habe den Code bisschen aufgeräumt, aber das problem besteht weiterhin.

Warum erhalte ich mit tfestimate nicht die erwartete Darstellung, die mit dem Alternativen Weg möglich ist.

Mein Verständnis ist, dass tfestimate, auch intern fft und die übertragungsfunction bildet

Bitte um Hilfe.


Code:


%% TFESTIMETE STARTING
clc; clear; close all

analysis = 'YAWRSWA';%
load('YAWRSWA.mat')    % YAW Rate
data       = YAWR;

fs         = FS;
WindowPara = [0.1881, -0.36923, 0.28702, -0.13077, 0.02488]; % Literatur
Resolution = 0.35;                                           % Literatur
Overlap    = 0.785;                                          % Literatur

NT = length(data)* (1/ fs);         % Länge des Zeitfensters in Sekunden
Abstand = 1/NT;                     % Abstand der einzelen Linien Resolution
N = ((0.5 * fs)/ Resolution)* 2;    % Blocklänge
Res = ((0.5 * fs)*2)/ length(data); % Abstand der einzelen Linien Resolution (max. possible)
% Do some argument validation
if N > length(data)
    disp([num2str(Res), ' Hz (max. Frequenzauflösung)'])
    error(message('Frequenzauflösung passt nicht zu Blocklänge'));
end
WindowSize = N;
nfft             = 2^nextpow2(length(data));
Overlap      = round(Overlap* WindowSize);
z                = 2*pi/WindowSize* (1:WindowSize-1);
Window      = WindowPara(1) + WindowPara(2)*cos(z)...
    + WindowPara(3)* cos(2*z) + WindowPara(4)*cos(3*z) + ...
    WindowPara(5)* cos(4*z);

 [magn, freq]   = tfestimate(SWA, YAWR, Window, Overlap, [], fs);
 [coh, freqcoh] = mscohere(SWA, YAWR, Window, Overlap, [], fs);
       
maxlim = 4;
figure;
subplot(3,1,1); hold all; grid on
p1 = plot(freq, real(magn));
xlabel('Frequenz [Hz]'); ylabel('Amplitude')
xlim([0 maxlim]); xticks(0:1:maxlim)
       
subplot(3,1,2); hold all; grid on
p1 = plot(freq, rad2deg(angle(magn)));
xlabel('Frequenz [Hz]'); ylabel('Phase [deg]')
xlim([0 maxlim]); xticks(0:1:maxlim)
       
 subplot(3,1,3); hold all; grid on
p1 = plot(freqcoh, coh);
xlabel('Frequenz [Hz]'); ylabel('Coherence [-]')
xlim([0 maxlim]); xticks(0:1:maxlim)
%  TFESTIMATE ENDING
       
%% ALTERNATIVE FFT and TRANSFER FUCNTION MANUALLY
x = SWA;
y = YAWR;
window = 1;
       
u =x1;
y =y1;
       
dt = 1/fs;
       
N = length(u);
U = fft(u);
Y = fft(y);
H1 = U .* Y ./ U.^2; % the H1 transfer function estimate
MagData = (abs(H1));
PhaseData = rad2deg(angle(H1));
       
N2 = floor(N/2);
FreqData = [0:N-1]/(N*dt); % Hz
       
figure;
subplot(211); plot(FreqData(1:N2),MagData(1:N2));
hold on; grid on
subplot(212); plot(FreqData(1:N2),PhaseData(1:N2));
hold on; grid on
       
n = length(x1);
c = fft(x1)/n;
amp = 2*abs(c);
amp(1) = amp(1)/2;
m = floor(N2/2);

% Ermittelt die letzte Frequenz mit mind. 10% Maximalamplitude!
frequenz=linspace(0,(m-1)*fs/n,m);
f_end=max(frequenz);
        for i=1:m
            i=m+1-i;
            if amp(i)>0.1*max(amp)
                f_end=frequenz(i);
                break
            end
        end
       
        for i=1:N2
            if FreqData(i)>f_end
                MagData=MagData(1:i);
                PhaseData=PhaseData(1:i);
                FreqData=FreqData(1:i);
                N2=i;
                break
            end
        end
       
        low_lim=0.2; %Nicht größer als 0.2Hz setzen! Sonst fehlen natürlich die Gain_0_2Hz Werte!
        for p2=1:N2-1
            if FreqData(p2)<=low_lim && FreqData(p2+1)>low_lim
                break
            end
        end
       
MagData    = MagData(p2(1):N2);
FreqData   = FreqData(p2(1):N2);
PhaseData  = PhaseData(p2(1):N2);
N2         = N2-p2(1);
figure
subplot(1,2,1); hold all; grid on
h = plot(FreqData(1:N2),MagData(1:N2),'linewidth',2); hold all
axis([low_lim maxlim -inf inf]);
smooingdata =smoothdata(MagData, 'movmedian');
h = plot(FreqData(1:N2),smooingdata(1:N2),'linewidth',2);
xlabel('Frequenz [Hz]'); ylabel('Magnitude [deg]')
subplot(1,2,2); hold all; grid on
h=plot(FreqData(1:N2),PhaseData(1:N2),'linewidth',2); hold all
axis([low_lim maxlim -inf inf]);
smooingdata =smoothdata(PhaseData, 'movmedian');
h=plot(FreqData(1:N2),smooingdata(1:N2),'linewidth',2);
xlabel('Frequenz [Hz]'); ylabel('Phase [deg]')
 
Private Nachricht senden Benutzer-Profile anzeigen
 
ThKo26
Themenstarter

 Forum-Century

Forum-Century

Beiträge: 184
Anmeldedatum: 21.09.17
Wohnort: ---
Version: 2015b
     Beitrag Verfasst am: 17.06.2020, 11:20     Titel:
  Antworten mit Zitat      
Hallo zusammen,

die Frage mit der Phase konnte ich klären, ich habe nochmal die Rohdaten genauer angeschaut und es verläuft gegenphasig.

Weiterhin ist die Frage offe, wieso die Ergebnisse zwischen den beiden Vorgehensweisen unterschiedlich sind.

Danke

Viele Grüße
Private Nachricht senden Benutzer-Profile anzeigen
 
ThKo26
Themenstarter

 Forum-Century

Forum-Century

Beiträge: 184
Anmeldedatum: 21.09.17
Wohnort: ---
Version: 2015b
     Beitrag Verfasst am: 18.06.2020, 13:28     Titel:
  Antworten mit Zitat      
hallo zusammen,

die antwortet lautet:

Code:
maxlim = 4;
figure;
subplot(3,1,1); hold all; grid on
p1 = plot(freq, abs(magn)); % --> not real!! must to be abs!!!
xlabel('Frequenz [Hz]'); ylabel('Amplitude')
xlim([0 maxlim]); xticks(0:1:maxlim)


Viele Grüße
Private Nachricht senden Benutzer-Profile anzeigen
 
Neues Thema eröffnen Neue Antwort erstellen



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 | goMatlab RSS Button 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.