SToNCalculator.m

Calculate your signal to noise ratio

function out = SToNCalculator(signal, PlotName)

    plot1Title = 'Signal, Fit, Noise, and Mean Value';
    plot2Title = 'Spectrum';

    if(nargin > 1)
        plot1Title = [PlotName ' ' plot1Title];
        plot2Title = [PlotName ' ' plot2Title];
    end

    %k = fir1(500, 0.1);
    %signal = filtfilt(k, 1, signal);

    fitFunction = @(x, xdata) exp( -1 * xdata ./ x(4)) .* x(1) .* cos( 2 * pi * 20 * xdata + x(2) ) + x(3) + x(5) .* xdata;

    time = (0:length(signal)-1)/1000;

    fitOptions = optimset('TolFun',1E-30,'TolX',1E-10,'MaxFunEvals',1E4,'MaxIter',1E3);
    fitValues = lsqcurvefit( fitFunction, [1 0 0 1000 0], time, signal, [0 -pi -inf 0 -1], [inf pi inf inf 1], fitOptions );

    noise = signal - fitFunction(fitValues, time);

    figure
    plot(time, signal, 'linewidth', 3)
    hold on
    plot(time, fitFunction(fitValues, time),'r');
    plot(time, noise + mean(signal), 'g');
    plot([time(1) time(end)], [mean(signal) mean(signal)]); 
    title(plot1Title)
    xlabel('Time (s)');
    ylabel('Signal (V)'); % 2013/11/21 was: 'Y(t) (V)'
    legend('Signal', 'Fit', 'Noise', 'Mean');

    figure 
    [Psignal, w] = pwelch(signal);
    [Pnoise, w] = pwelch(noise);
    w = w * 500 / max(w);
    frequencyAxis = w * 500 / max(w);
    semilogy(frequencyAxis, Psignal,'b');
    hold on;
    semilogy(frequencyAxis, Pnoise, 'g');
    xlabel('Frequency (Hz)');
    ylabel('Magnitude'); % 2013/11/21 was: 'Magnitude (dB)'
    title(plot2Title);

    out = {};
    out.Amplitude = fitValues(1);
    out.RMSAmplitude = out.Amplitude / sqrt(2);
    out.Phase = fitValues(2);
    out.Offset = fitValues(3);
    out.BleachingTimeConstant = fitValues(4);
    out.Tilt = fitValues(5);
    out.Signal = signal;
    out.FitSignal = fitFunction(fitValues, time);
    out.Time = time;
    out.Noise = noise;
    out.FrequencyAxis = frequencyAxis;
    out.SignalMagnitudeSpectrum = Psignal;
    out.NoiseMagnitudeSpectrum = Pnoise;
    out.RMSE = sqrt(mean(noise.^2));  % 2014/11/20 was: (noise/fitValues(1))
    out.AdjustedRMSE = out.RMSE / exp(-60*10 / out.BleachingTimeConstant);
    out.SNR = 20 * log(out.RMSAmplitude / out.RMSE) / log(10);
    out.AdjustedSNR = 20 * log(out.RMSAmplitude / out.AdjustedRMSE) / log(10);

end

% by SCW, Spring 2010