%% 02/03/2012
% Prove per curve di potenza - resistore 0-100 kohm
% prove a +-0.5g
% R=10,20,30,40,50,60,70,80,90,97 kohm
% m=0,2.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=[1 2 3 4 5 6 7 8 9 9.7]*1e4; %ohm

% vettore masse
m=[0 2.6]; %g


%lista  file di input e rispettive funzioni del datasheet
lista=struct('file',{[21:30];[31:40]},...
             'funzione',{'-0.0446*x.^2+0.2897*x';'-0.0353*x.^2+0.2687*x'});

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



%% Ricerca delle frequenze proprie

% con trasformata (ris in freq df=20Hz!)
f=[];
for i=[21 31]
    dati=importdata(['dati/scope_' int2str(i) '.csv'],',',2);
    acce=dati.data(:,3);
    
    N=length(acce);
    dt=(-dati.data(1,1)+dati.data(2,1));
    fsamp=1/dt;
    df=fsamp/N;
    
    %trasformata
    TF=fft(acce)./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=[21 31]
    dati=importdata(['dati/scope_' int2str(i) '.csv'],',',2);
    acce=dati.data(:,3);
    dt=(-dati.data(1,1)+dati.data(2,1));
    fsamp=1/dt;
        
    [ac tau]=xcorr(acce,'unbiased');
%     figure %plot dell'autocorrelazione
%     plot(tau,ac,'r')
%     xlabel('tau   [s]')
%     ylabel('autocorrelazione   [EU^2]')
        
    [picchi,val]=pickpeak(ac(1:(end+1)/2),2,50);
    periodo=abs(mean(diff(picchi)/fsamp));
    freq=1/periodo;
    
    f=[f freq];
end

% frequenze proprie per ciascuna delle 4 configurazioni
f=[120.4 75.0]

x=[0:.1:30];
  
%% plot curve datasheet e risultati sperimentali
for i=1:length(lista)
    
    clear dV
    clear P
    clear dati
    for k=1:length(lista(i).file)
        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),2);
%         dV(k)=mean(val); % PEAK VALUE
        [picchi_ac,val_ac]=pickpeak(acc-mean(acc),2,50);
        [picchi_amp,val_amp]=pickpeak(amp-mean(amp),2,50);

        dV(k)=sqrt(mean(harv)^2+std(harv)^2); %RMS
        
        ACC(i,k)=mean(val_ac);
        AMP(i,k)=mean(val_amp);
        
        P(k)=dV(k).^2./R(k)*1000;
        
%         figure(10+i) % output debug
%         plot(harv-mean(harv))
%         hold on
        
    end
    

    f0=f(i); %frequenza fondamentale
        
    %funzione curva datasheet
    funz=inline(lista(i).funzione);
    figure(1)    
    subplot(1,2,i)
    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),f0);
    title(titolo);
    xlabel('V [V_R_M_S] ')
    ylabel('P [mW]')
    
    
    figure(2)    
    hold on
    plot(dV,P,[colore(i) 'o-'])
    ylim([0,inf])
    grid on
    title('potenze sperimentali')
    xlabel('V [V_R_M_S]')
    ylabel('P [mW]')
    
    figure(3)
    hold on
    plot(R(1:length(lista(i).file)),ACC(i,1:length(lista(i).file)),[colore(i) 'o-'])
    grid on
    title('accelerazioni')
    xlabel('R [\Omega]')
    ylabel('a [m/s^2]')
    
    figure(4)
    hold on
    plot(R(1:length(lista(i).file)),AMP(i,1:length(lista(i).file)),[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;
    
%     dV  % output debug
%     pause 
%     tens(i,1:length(lista(i).file))=dV;
%     pot(i,1:length(lista(i).file))=P;
end


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

% % salvataggio per cfr varie prove
% save tensioni4.dat tens -ascii
% save potenze4.dat pot -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