function [Amplitude,Frequency,DC] = noisin(data,Fs,method)
% noisin.m - Estimate the amplitude, frequency, and DC offset of a sampled noisy sinusoid.
%
% Syntax: [Amplitude,Frequency,DC] = noisin(data vector,sampling rate,'method')
%
% The frequency output will be the frequency component with the greatest amplitude (Not including f=0).
%
% DC amplitude is evaluated from the signal Fourier transform at f=0.
%
% The amplitude estimate is evaluated in one of two ways: The RMS method is better
% with low sampling rates or very high S/N, the amplitude spectrum method is better with lower S/N
% or non-random noise. The spectrum method is generally more sensitive to the frequency estimate.
% While the spectrum method is more accurate in most situations, the sampling rate/frequency ratio
% and the number of waveform cycles in the data set are positively correlated to the accuracy of
% the results -- the spectral method of amplitude estimation being more sensitive to both factors
% than the RMS method. The program can become significantly unreliable if the sinusoid is composed
% of two or more components with almost equal amplitudes and/or frequencies. In this regard, some
% prior knowledge about the physical process that produced the real-valued sinusoid, or about its
% frequency structure, may be indispensable.
%
% One way to decide on the amplitude estimation method is to generate an artificial noisy sinusoid
% of known dominant frequency and amplitude, and with roughly the same amplitude, frequency,
% number of samples, sampling rate, and noise characteristics as the sinusoid in question, and then
% to compare the accuracy of the results of the RMS and of the spectrum methods, obtained for the
% artificial sinusoid.
%
% Set 'method' = 'rms' to select the RMS method, or 'as' to select the amplitude spectrum method.
%
%
% For a very noisy signal, or for non-random noise, try filtering around the frequency
% component of interest, and compensating for any attenuation in the passband before passing
% the filtered sinusoid to noisin.m, in order to improve the amplitude estimate. The frequency
% estimate is of course contingent on the frequency satisfying the Nyquist criterion, i.e.,
% that the sampling rate be at least double the frequency component to be estimated.
%
% Important Note 1: Do not pad the raw data with zeros. Include as much of
% the continuous steady-state portion of the waveform as possible, and as
% little of the transient or flat portions as possible. Truncate if necessary,
% but avoid concatenating waveforms.
%
% Important Note 2: This program will give erroneous results if the main frequency
% is modulated, or multiplied, by another. A sinusoidal waveform is a waveform
% formed from the superposition of several sine terms, not from their multiplication.
%
% See the code for more details.
%
% Author: Hazem Saliba Baqaen (Hazem@brown.edu)
% Thanks to Sami Deeb and Eric LePage for suggestions and feedback.
% Designed on MATLAB 6.5.
[r,c] = size(data); % Check dimensions of data set
if (r~=1) & (c~=1)
error('Data must be a one-dimensional vector')
elseif (r~=1)
data = data';
end
clear r,c;
if Fs ~= abs(Fs) % Make sure Fs is real and positive
error('The sampling rate must be a real positive number')
end
[DC,Frequency,MaxAmp,datanodc] = ampspectrum(data,Fs);
if strcmpi('rms', lower(method))
% Amplitude of a noisy sinusoid (no DC) approximately= (RMS value of noisy sinusoid / RMS value of A = 1 noiseless sine wave)
% This method of estimating the amplitude relies only on the time-domain signal representation and is independent of the
% frequency of the sinusoid (Except in that the frequency is needed to make sure only an integer number of cycles are used
% in the RMS amplitude estimate).
% New 'datanodc' data set (with the estimated DC component subtracted out) comes from ampspectrum subfunction
Periodn = Fs/Frequency; % Period in samples
datanodc = datanodc(1:fix(fix(length(datanodc)/Periodn)*Periodn)); % Make sure only an integer number of cycles is used in the RMS estimate
Amplitude = (norm(datanodc)/sqrt(length(datanodc)))/sqrt(0.5);
elseif strcmpi('as', lower(method))
% Alternative (Spectral) method of estimating the amplitude (using frequency-domain amplitude estimate).
Amplitude = MaxAmp;
else
Amplitude = MaxAmp; % The spectrum method of amplitude estimation is the default method
end
function [DC,Frequency,MaxAmp,datanodc] = ampspectrum(data,Fs)
% Obtain simple amplitude spectrum estimate, then find DC amplitude, subtract
% from signal, then find non-DC main frequency and its amplitude.
% Adapted from: http://www.mathworks.com/support/tech-notes/1700/1702.html
% Define a raised-cosine-squared windowing of data for increasing smoothness and
% improving estimates resulting from the fft of the data
raisedcosinesqwindow = (2/3)*(1-cos(2*pi*[0:length(data)-1]/length(data))).^2;
% Use next highest power of 2 greater than or equal to length(x=data) to calculate FFT.
%% (To take advantage of speeded fft algorithm. Pads with zeros so that length(data) is a radix-2 number)
NFFT= 2^(nextpow2(length(data))); %(length(data));
% Take the Fourier Transform of the waveform
FFTX = fft(data.*raisedcosinesqwindow,NFFT);
% % Calculate the number of unique points
% NumUniquePts = ceil((NFFT+1)/2);
%
% % FFT is symmetric, throw away second half
% FFTX = FFTX(1:NumUniquePts);
%
% % Take the magnitude of fft of x
% MX = abs(FFTX);
%
% % Scale the fft so that it is not a function of the
% % length of x
% MX = MX/length(data);
%
% % Multiply by 2 because you
% % threw out second half of FFTX above
% MX = 2*MX;
%
% % DC Component should be unique.
% MX(1) = MX(1)/2;
%
% % Nyquist component should also be unique.
% if ~rem(NFFT,2)
% % Here NFFT is even; therefore, Nyquist point is included.
% MX(end) = MX(end)/2;
% end
%
% % This is an evenly spaced frequency vector with
% % NumUniquePts points.
% f = (0:NumUniquePts-1)*Fs/NFFT;
% figure
% plot(f,MX);
% title('Amplitude Spectrum of Signal');
% xlabel('Frequency ');
% ylabel('Amplitude');
% grid
% Find DC amplitude, subtract from signal, then find non-DC main frequency.
% Frequency estimate is improved if taken from the amplitude spectrum of
% the autocorrelation function of the signal data, rather than from the
% amplitude spectrum of the signal itself.
% New: The frequency and the amplitude estimates can be improved further
% still by using the eigenvectors of the autocorrelation matrix method.
% See rooteig.m and rootmusic.m
DC = sign(real(FFTX(1)))*abs(FFTX(1))/length(data);
data = data - DC;
datanodc = data; % Save copy of data with DC component subtracted out under another name
acorrdata = xcorr(data); % Autocorrelation function for the signal (Without DC)
acorrdata = [acorrdata zeros(size(acorrdata)) zeros(size(acorrdata))]; % Triple the length padding with zeros to improve frequency estimate
NFFT_acorr = length(acorrdata);
FFTX_acorr = fft(acorrdata);
NumUniquePts_acorr = ceil((NFFT_acorr+1)/2);
FFTX_acorr = FFTX_acorr(1:NumUniquePts_acorr);
MX_acorr = abs(FFTX_acorr);
MX_acorr = MX_acorr/length(acorrdata);
MX_acorr = 2*MX_acorr;
MX_acorr(1) = MX_acorr(1)/2;
if ~rem(NFFT_acorr,2)
MX_acorr(end) = MX_acorr(end)/2;
end
f_acorr = (0:NumUniquePts_acorr-1)*Fs/NFFT_acorr;
% Find frequency component with maximum amplitude (Excluding DC component)
[MaxAmpAcorr, MaxAmpAcorrindex] = max(MX_acorr);
Frequency = f_acorr(MaxAmpAcorrindex); % This gives a close estimate (usually +/- 0.5 Hz) that allows us to put an upper cutoff for the more accurate frequency estimate (see below)
FirstEstimate = Frequency; % Remember first estimate of frequency
% % Trim spectrum around first frequency estimate, then use to get a closer
% % estimate using rooteig or rootmusic (currently not wroking properly)
% NFFT = length(data);
% FFTX = fft(data);
% NumUniquePts = ceil((NFFT+1)/2);
% FFTX = FFTX(1:NumUniquePts);
% f = (0:NumUniquePts-1)*Fs/NFFT;
% narrowbandindices = find(abs(f-Frequency)<2);
% %filtereddata = real(ifft([zeros(1,min(narrowbandindices)-1) FFTX(narrowbandindices) zeros(1,NFFT-max(narrowbandindices))],NFFT));
% filtereddata = real(ifft(FFTX(narrowbandindices),NFFT));
% Very accurate autocorrelation eigenvector estimation method: -------\
if length(data) <= 100000 % Use this estimation method only if the total number of samples is less than some max. This is to avoid OUT OF MEMORY errors and long waits. The max depends on p, among other factors.
p=40; % Twice the number of frequency components to estimate. Has to be large enough and even (see rooteig.m and rootmusic.m). This method is especially suited for line spectra.
[freqcompts,P] = rooteig(data,p,Fs);
freqcompts = freqcompts(freqcompts >= 0); % We are dealing with a real sinusoid (not complex); drop negative frequencies
freqcompts = freqcompts(1:length(P)); % Truncate freqcompts to match number of elements in power vector
P = P(freqcompts <= (Frequency + 15.0)); % Set upper limit on frequency, find associated power values
freqcompts = freqcompts(freqcompts <= (Frequency + 15.0));
P(freqcompts == 0) = 0; % The power (and amplitude) of the DC component should be zero because it has already been subtracted out.
%% Force it to zero (without modifying the signal) in order to prevent it from masking any small sinusoidal components.
[maxP,maxPindex] = max(P);
Frequency = freqcompts(maxPindex); % Dominant frequency estimate
MaxAmp = sqrt(2*maxP); % Amplitude of dominant frequency
% stem(freqcompts,sqrt(2*P),'+');grid
end
% --------------------------------------------------------------------/
% Revert to first estimates if the eigenvector method fails
if (isempty(Frequency) || (abs(FirstEstimate - Frequency) > 0.5)) | (length(data) > 100000)
% Older method of amplitude estimation
%
% Find the amplitude of the main frequency (Other than at f=0)
data = [data.*raisedcosinesqwindow zeros(size(data))]; % Raised-cosine-squared 'windowing' of limited-time signal to reduce artifacts and improve estimate...
%...and then padding with zeros to increase the size of NFFT.
NFFT = length(data);
FFTX = fft(data);
NumUniquePts = ceil((NFFT+1)/2);
FFTX = FFTX(1:NumUniquePts);
MX = abs(FFTX);
MX = MX/length(data);
MX = 2*MX;
MX(1) = MX(1)/2;
if ~rem(NFFT,2)
MX(end) = MX(end)/2;
end
f = (0:NumUniquePts-1)*Fs/NFFT;
MaxAmp = 2*max(MX);
Frequency = FirstEstimate;
end
