%% 17/02/2012
% Prove per curve di potenza - seconda parte
% prove a +-1g e +-0.5g
% R=1,2,3,4,5,6,7,8,9,9.86 kohm
% m=0,2.4,7.8,15.6 grammi

close all
clear all
clc

%% Struttura dei file di input:
% tempo [s]
% channel 2: spostamento tip laser [V]
% channel 3: accelerometro su shaker [V]
% channel 4: harvester [V]

%sensibilità 
laser_sens= 2; %V/mm
acc_sens=100e-3/9.81; %V/(m/s^2)

% vettore resistenze
R=[9.86 9 8 7 6 5 4 3 2 1]*1e3; %ohm

% vettore masse
m=[0 2.4 7.8 15.6]; %g


%lista  file di input e rispettive funzioni del datasheet
lista=struct('file',{[2:11];[12:21];[23:32];[33:42]; 
                     [44:53];[54:63];[66:75];[76:85]},...
             'funzione',{'-0.0480*x.^2+0.5338*x';'-0.0446*x.^2+0.2897*x';...
                         '-0.0353*x.^2+0.5861*x';'-0.0353*x.^2+0.2687*x';...
                         '-0.0412*x.^2+1.0892*x';'-0.0334*x.^2+0.5256*x';...
                         '-0.0413*x.^2+1.2357*x';'-0.0408*x.^2+0.6919*x'});

 
%figure('position',[100,100,800,600])
% vettore colori
col=['brbrbrbrbrbrbrbr'];



%% Ricerca delle frequenze proprie

% con trasformata (ris in freq df=20Hz!)
f=[];
for i=[1 22 43 65]
    dati=importdata(['dati/scope_' int2str(i) '.csv'],',',2);
    spost=dati.data(:,2);
    
    N=length(spost);
    dt=(-dati.data(1,1)+dati.data(2,1));
    fsamp=1/dt;
    df=fsamp/N;
    
    %trasformata
    TF=fft(spost)./N;
    modulo=abs(TF(1:N/2+1))*2;
    modulo(1)=abs(TF(1));
    modulo(N/2+1)=abs(TF(N/2+1));
    freq=0:df:fsamp/2;

    %plot(freq,modulo)
    
    [massimo,pos]=max(modulo);
    f=[f freq(pos)];
end


% con autocorrelazione
f=[];
for i=[1 22 43 65]
    dati=importdata(['dati/scope_' int2str(i) '.csv'],',',2);
    spost=dati.data(:,2);
    dt=(-dati.data(1,1)+dati.data(2,1));
    fsamp=1/dt;
        
    [ac tau]=xcorr(spost,'unbiased');
%     figure %plot dell'autocorrelazione
%     plot(tau,ac,'r')
%     xlabel('tau   [s]')
%     ylabel('autocorrelazione   [EU^2]')
        
    [picchi,val]=pickpeak(ac(1:(end+1)/2),3);
    periodo=abs(mean(diff(picchi)/fsamp));
    freq=1/periodo;
    
    f=[f freq];
end

% frequenze proprie per ciascuna delle 4 configurazioni
f
x=[0:.1:30];
  
%% plot curve datasheet e risultati sperimentali
for i=1:length(lista)
    
    clear dV
    clear dati
    for k=1:length(R)
        dati=importdata(['dati/scope_' int2str(lista(i).file(k)) '.csv'],',',2);
        harv=dati.data(:,4);
        acc=dati.data(:,3)./acc_sens;
        amp=dati.data(:,2)./laser_sens;
%         [picchi,val]=pickpeak(harv-mean(harv),3,50);
%         dV(k)=mean(val);  %%PEAK VALUE
        [picchi_ac,val_ac]=pickpeak(acc-mean(acc),3,50);
        [picchi_amp,val_amp]=pickpeak(amp-mean(amp),3,50);

        ACC(i,k)=mean(val_ac);
        AMP(i,k)=mean(val_amp);

        dV(k)=sqrt(mean(harv)^2+std(harv)^2); %RMS  
        
    end
    
    
    P=dV.^2./R*1000; %mW

    f0=f((i+mod(i,2))/2); %frequenza fondamentale
        
    %funzione curva datasheet
    funz=inline(lista(i).funzione);
    figure(1)    
    subplot(2,2,(i+mod(i,2))/2)
    hold on
    plot(x,funz(x),col(i),'linewidth',2)
    plot(dV,P,[col(i) 'o'])
    ylim([0,inf])
    grid on
    titolo=sprintf('m= %.1f g, f_0= %.1f Hz',m((i+mod(i,2))/2),f0);
    title(titolo);
   % xlabel('V [V]')
    xlabel('V [V_R_M_S]')
    ylabel('P [mW]')
    
    
    colore=['ggyyrrbb'];
    
    figure(2)    
    hold on
    plot(dV,P,[colore(i) 'o-'])
    ylim([0,inf])
    grid on
    title('potenze sperimentali')
   % xlabel('V [V]')
    xlabel('V [V_R_M_S]')
    ylabel('P [mW]')
    
    figure(3)
    hold on
    plot(R,ACC(i,:),[colore(i) 'o-'])
    grid on
    title('accelerazioni')
    xlabel('R [\Omega]')
    ylabel('a [m/s^2]')
    
    figure(4)
    hold on
    plot(R,AMP(i,:),[colore(i) 'o-'])
    grid on
    title('ampiezze')
    xlabel('R [\Omega]')
    ylabel('a [mm]')
    
    risposta(i,2)=AMP(i,1)./ACC(i,1);
    risposta(i,1)=f0;

%     tens(i,:)=dV;
%     pot(i,:)=P;
    
end

% % salvataggio per cfr varie prove
% save tensioni2.dat tens -ascii
% save potenze2.dat pot -ascii




% output risultati per cfr FRF
save risposta.txt risposta -ascii

% salvataggio EPS
figure(1)
print -painters -depsc img/curve_potenze.eps
print -painters -dpng img/curve_potenze.png
figure(2)
print -painters -depsc img/potenze.eps
print -painters -dpng img/potenze.png
figure(3)
print -painters -depsc img/accelerazioni.eps
print -painters -dpng img/accelerazioni.png
figure(4)
print -painters -depsc img/ampiezze.eps
print -painters -dpng img/ampiezze.png