function [mag, mag_dB, fv, phase_rad, phase_grad] = FFT_bode( signal, nfft, fa) ; 
 
% Input: 
% Signal im Zeitbereich 
% nfft = Anzahl Messwerte für fft 
% wenn nfft > length(sig) -> fft(sig,nfft) führt Zeropadding durch 
% fa = Abtastfreq. 
% Output: 
% Magnitude des Spektrums linear und dB skaliert 
% Frequenzvektor fv in [Hz] von 0...fa/2 
fa = 4000;
d = evalin('base','Y_eff')%load('vc50_vf18_ae4_f3_1__s3.txt');%)%
signal = d;
nfft = length(signal);
% un-,gerade Anzahl Messwerte? 
if mod(nfft,2) == 0; 
    k = (nfft/2) + 1; 
else 
    nfft = nfft + 1; 
    k = (nfft/2) + 1; 
end 

fn = fa/2; % Nyquistfreq. 
df = fa/nfft; % Frequenzauflösung des Spektrums 
% Frequenzvektor: Darstellung bis Nyquistfreq. 
fv = 0: df : fn; 

sig = signal(:); 

% Signal transformieren 
Fy = fft(sig,nfft); 
  
%   Betrag - nur positives Freq.spektrum   
Fy_pos0 = abs(Fy(1:k)); 
  
%   Skalierung   
mag = [Fy_pos0(1)/nfft ;Fy_pos0(2:k-1)/(nfft/2);Fy_pos0(k)/nfft]; 

% Umrechnung in dB 
mag_dB = 20*log10(mag + eps); % eps = kleine Konstante zur Vermeidung von log(0) 

% Berechung der Phase in rad/s 
phase_rad = angle(Fy(1:k)); 

% Umrechnung in grad 
phase_grad = phase_rad * (180/pi) ;


subplot(2,1,1);
semilogx(fv,mag_dB) 
xlabel('Freuqenz in Hz'); 
ylabel('Magnitude in dB');
xlim([min(fv) max(fv)])
grid on;

subplot(2,1,2); 
semilogx(fv,phase_grad); 
ylabel('Phase [°]'); 
xlabel('Frequenz [rad/s]');
xlim([min(fv) max(fv)]); 
grid on;
end