function CpRatio = CpRatio_16Aug2017(lightWaveform) % Script for calculating Cp Ratio: Cp(DUT)/Cp(Square wave) % Function arguement % lightWaveform: A two-column array of light output and time in equal increments % Column 1: time in units of seconds % Column 2: relative light output time = lightWaveform(:,1); deltaTime = time(2)-time(1); Light = lightWaveform(:,2); FS = 1/deltaTime percentFlicker = 100*(max(Light)-min(Light))/(max(Light)+min(Light)); DC = mean(Light); % Determine frequency of maximum spectral component in light waveform Lightf = fft(Light)/length(Light); if (rem(length(time),2)==0) % even length %freq = ifftshift(linspace(-1/(2*deltaTime),1/(2*deltaTime)-1/(time(end)-time(1)),length(Lightf))); freq = ifftshift((0:(length(time)-1))*(FS/length(time))-FS/2); % Equivalent to line above, but not roundoff error else % odd length freq = ifftshift(linspace(-1/(2*deltaTime),1/(2*deltaTime),length(Lightf))); end Lightf = filtfilt([.2,.2,.2,.2,.2],1,Lightf); % Filter as a way of integrating closely spaced freq components [Max,index] = max(abs(Lightf(6:floor(length(Light)/2)-1))); % Skip the first 5 elements as they contain dc leakage freqOfMaximum = freq(index+5) A = DC*percentFlicker/100; dutyCycle = 50; LightSqWave = (A*square(2*pi*freqOfMaximum*time,dutyCycle))+DC; Window = window(@hann,length(Light)); % To minimize finite sample length effects Light = Light.*Window; LightSqWave = LightSqWave.*Window; A_DUT = sum(Light.*exp(1i*freqOfMaximum*2*pi*time)); % Single component DFT A_sqWave = sum(LightSqWave.*exp(1i*freqOfMaximum*2*pi*time)); % % Single component DFT CpRatio = abs(A_DUT)/abs(A_sqWave);