close all
load Input_data
size = 12;
N = 4096;  
Fs = 500;   
Nyquist = 250;   

freq_axis = linspace(0,Fs,N);     %N points between 0 and Fs

x = detrend(alfa(1:N,2));
xx = detrend(Input(1:N,2));
h = spectrum.welch;                    % Create a Welch spectral estimator. 
[Px, w_welch] = pwelch(x,1024,32,1024,Fs);      % Calculate the PSD 
[Pxx, w_welch] = pwelch(xx,1024,32,1024',Fs);   % Calculate the PSD 
semilogy(w_welch, Px./Pxx, 'r','LineWidth',1)

hold on
y = abs(fft(x)).^2;
yy = abs(fft(xx)).^2;
plot(freq_axis,y./yy,'b-')
xlim([0 50])
title('Power Spectral Density','FontSize', size, 'VerticalAlignment', 'bottom')
set(gca, 'Fontsize', size)

ylabel('dB/Hz','FontSize',size)
xlabel('Hz','FontSize',size, 'VerticalAlignment', 'middle')

A = 1.65;
a1 = 55;
a2 = 605;
B = 32;
b1 = 27.5;
b2 = 55;
c1 = 32;
c2 = 3;
Ve = 1.80;  %operating point
Vi = 2.93;  %operating point

%derrivative
%f' = -fmax*(exp()*(1/sigma))/(1+exp())^2;
Ve_exp = exp((Ve-7)/-2);
Vi_exp = exp((Vi-7)/-2);

q1 = 0.5*50*Ve_exp/(1+Ve_exp)^2;
q2 = 0.5*50*Vi_exp/(1+Vi_exp)^2;

% s = sym('s');                %creates symbolic variable
% 
% num = A*(a2-a1)*(b1+s)*(b2+s);        %numerator
% den = (a1+s)*(a2+s)*(b1+s)*(b2+s);    %denominator
% 
% 
% N = sym2poly(num);         %coefficients of a polynomial
% D = sym2poly(den);

N = [907.5  74868.75   1372593.75];
D = [1   742.5    89237.5   3743437.5    50328437.5];

K= c1*c2*q1*q2*A*(a2-a1)*B*(b2-b1);
D(5) = D(5) + K;      

w = linspace(0,314);    %frequency axis in to radians
f = w/(2*pi);           %convertion to Hertz
h = freqs(N,D,w);       %frequency response in s-domain 
mag = abs(h);           %magnitude of frequency response 
plot(f,mag.^2,'k','LineWidth',1.5);

legend('Welch','|fft(V)|^2', 'analytic')







