function [Z,fr,varA, varPhi,stats,Wz,c_u,c_i,ts]=fc_impedance_cwtft(voltage,current,f_min,f_max,fs,decimate_factor,f0,nb)
    %f0=5;
    if (decimate_factor > 1)
        u_fc = decimate(voltage,decimate_factor);
        i_fc = decimate(current,decimate_factor);
        fs = fs/decimate_factor;
    else
        u_fc = voltage;
        i_fc = current;
    end
    
    
    %f0=1;
    MorletFourierFactor = 4*pi/(f0+sqrt(2+f0^2));

    WAV={'morl',f0};
    %nb = 23;

    s0 = fs/(f_max * MorletFourierFactor);
    mx = fs/(f_min * MorletFourierFactor);
    ds = log(mx/s0)/log(2)/(nb-1);

    %s0 = (0.002*MorletFourierFactor);
    %smax = (3/100*MorletFourierFactor);
    %ds =0.2;
    %nb = ceil(log2(smax/s0)/ds+1);
    %scales = struct('s0',s0,'ds',0.5,'nb',nb,'type','pow','pow',2);
    scales = struct('s0',s0,'ds',ds,'nb',nb,'type','pow','pow',2);

    disp('Starting CWT..');
    c_u = cwtft(u_fc,'wavelet',WAV,'scales',scales);
    c_i = cwtft(i_fc,'wavelet',WAV,'scales',scales);
    disp('CWT Done..');
    Wz=c_u.cfs./c_i.cfs;
    
    
    %Cone of Influence for Morlet wavelet
    %http://paos.colorado.edu/research/wavelets/bams_79_01_0061.pdf
    %A Practical Guide to Wavelet Analysis, Christopher Torrence and Gilbert P. Compo
    s=c_u.scales;
    ts=ceil(sqrt(2)*s);
    %disp(['ts=' num2str(length(ts))]);
    varA = zeros(size(ts));
    varPhi = zeros(size(ts));
    Z = zeros(size(ts));
    disp(['Calculating stat for ' num2str(length(ts)) ' scales...']);
    for ii=1:length(ts)
        if (size(Wz,2) <= ts(ii))
            varA(ii) = NaN;
            varPhi(ii) = NaN;
            Z(ii) = NaN;
            fr(ii) = NaN;
            
            stats.angle(ii).median = NaN;
            stats.angle(ii).prctile10 = NaN;
            stats.angle(ii).prctile90 = NaN;
            
            stats.abs(ii).median = NaN;
            stats.abs(ii).prctile10 = NaN;
            stats.abs(ii).prctile90 = NaN;
            
            stats.real(ii).median = NaN;
            stats.real(ii).prctile10 = NaN;
            stats.real(ii).prctile90 = NaN;
            
            stats.imag(ii).median = NaN;
            stats.imag(ii).prctile10 = NaN;
            stats.imag(ii).prctile90 = NaN;
            
            %disp(ii);
        else
            varA(ii) = var(abs(Wz(ii,ts(ii):end-ts(ii))));
            varPhi(ii) = var(angle(Wz(ii,ts(ii):end-ts(ii))));
            Z(ii) = mean(real( Wz(ii,ts(ii):end-ts(ii)) ) ) ...
                + 1i*mean(imag( Wz(ii,ts(ii):end-ts(ii)) ) );
            
            %stats.angle(ii).median = median(angle(Wz(ii,ts(ii):end-ts(ii))));
            %stats.angle(ii).prctile10 = prctile(angle(Wz(ii,ts(ii):end-ts(ii))),10);
            %stats.angle(ii).prctile90 = prctile(angle(Wz(ii,ts(ii):end-ts(ii))),90);
            
            %stats.abs(ii).median = median(abs(Wz(ii,ts(ii):end-ts(ii))));
            %stats.abs(ii).prctile10 = prctile(abs(Wz(ii,ts(ii):end-ts(ii))),10);
            %stats.abs(ii).prctile90 = prctile(abs(Wz(ii,ts(ii):end-ts(ii))),90);
            
            stats.real(ii).median = median(real(Wz(ii,ts(ii):end-ts(ii))));
            stats.real(ii).prctile10 = prctile(real(Wz(ii,ts(ii):end-ts(ii))),10);
            stats.real(ii).prctile90 = prctile(real(Wz(ii,ts(ii):end-ts(ii))),90);
            
            stats.imag(ii).median = median(imag(Wz(ii,ts(ii):end-ts(ii))));
            stats.imag(ii).prctile10 = prctile(imag(Wz(ii,ts(ii):end-ts(ii))),10);
            stats.imag(ii).prctile90 = prctile(imag(Wz(ii,ts(ii):end-ts(ii))),90);
            
            
            %stats.current.real(ii).mean = mean(real(c_i.cfs(ii,ts(ii):end-ts(ii))));
            %stats.current.real(ii).var = var(real(c_i.cfs(ii,ts(ii):end-ts(ii))));
            %stats.current.imag(ii).mean = mean(imag(c_i.cfs(ii,ts(ii):end-ts(ii))));
            %stats.current.imag(ii).var = var(imag(c_i.cfs(ii,ts(ii):end-ts(ii))));
            
            %stats.voltage.real(ii).mean = mean(real(c_u.cfs(ii,ts(ii):end-ts(ii))));
            %stats.voltage.real(ii).var = var(real(c_u.cfs(ii,ts(ii):end-ts(ii))));
            %stats.voltage.imag(ii).mean = mean(imag(c_u.cfs(ii,ts(ii):end-ts(ii))));
            %stats.voltage.imag(ii).var = var(imag(c_u.cfs(ii,ts(ii):end-ts(ii))));
            
            
           % [stats.real(ii).hist_bin stats.real(ii).hist_x] = hist(real(Wz(ii,ts(ii):end-ts(ii))), linspace(0, 5, 2000));
           % [stats.imag(ii).hist_bin stats.imag(ii).hist_x] = hist(imag(Wz(ii,ts(ii):end-ts(ii))), linspace(-2, 0.5, 2000));
           % 
           % [stats.abs(ii).hist_bin stats.abs(ii).hist_x] = hist(abs(Wz(ii,ts(ii):end-ts(ii))), linspace(0, 5, 2000));
           % [stats.angle(ii).hist_bin stats.angle(ii).hist_x] = hist(angle(Wz(ii,ts(ii):end-ts(ii))), linspace(-2, 0.5, 2000));
            
           % [stats.real(ii).muhat stats.real(ii).sigmahat stats.real(ii).muci stats.real(ii).signaci] = normfit(real(Wz(ii,ts(ii):end-ts(ii))));
           % [stats.imag(ii).muhat stats.imag(ii).sigmahat stats.imag(ii).muci stats.imag(ii).signaci] = normfit(imag(Wz(ii,ts(ii):end-ts(ii))));
            
           % [stats.real(ii).tlocationscale] = mle(real(Wz(ii,ts(ii):end-ts(ii))), 'distribution', 'tlocationscale');
           % [stats.imag(ii).tlocationscale] = mle(imag(Wz(ii,ts(ii):end-ts(ii))), 'distribution', 'tlocationscale');
            %[stats.imag(ii).nakagami] = mle(imag(Wz(ii,ts(ii):end-ts(ii))), 'distribution', 'nakagami');
            
            fr(ii) = 1./(c_u.scales(ii)*MorletFourierFactor)*fs;
            if isnan(fr(ii))
                disp(c_u.scales(ii));
                disp(MorletFourierFactor);
                disp(fs);
            end
        end
    end
    
    %varA = var(abs(Wz),0,2);
    %varPhi = var(angle(Wz),0,2);
    
    %Z=mean(real(Wz),2)+1i*mean(imag(Wz),2);
    
    
    
    %fr = 1./(c_u.scales*MorletFourierFactor)*fs;
end
