SToNCalculator.m
From Course Wiki
Revision as of 18:44, 16 November 2015 by Josephinebagnall (Talk | contribs)
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