Hello, I am trying to analyze some data from an ADC output (although we don't have the ADC
yet and I'm using Matlab to generate a theoretical output).  We are putting a signal with a
frequency in and the output will contain the output, quantization noise, and harmonic spurs
(using a spectrum analyzer).  Everything is being viewed in the frequency domain with the
highest amplitude = 0 dB (using dBFS).  
I was having some trouble calculating the Signal to Noise Ratio (ratio of the original signal
to only the quantization noise).  The data that I have will contain the input frequency,
quantization noise, and harmonic spurs.  I need to get rid of the harmonic spurs and input
frequency.  There is some spectral leakage from the harmonic spurs and input frequency so I
used a matlab algorithm to find the spurs and add (in the linear scale, not logarithm) the bins
to the left and right of it to get the correct amplitude power.

I'm having trouble making an algorithm to get rid of the spurs and the input frequency spur.  I
can't just subtract the spurs away because I have to account for the quantization noise.  If I
could find a way to calculate the SNR then I believe I will be able to calculate the SINAD
(ratio of the signal to all the harmonics.  Would only need to get data with a subtracted input
frequency), THD (need a ratio of the input signal to the harmonic spurs), and SFDR (already
have this done, just subtracted the input frequency spur from the second spur).  
Here is what I have done so far:

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:1:L-1)*T;                % Time vector

fin = 100; %input frequency
amplitude = 10;  %amplitude of input signal

x = amplitude*sin(2*pi*fin*t); %input signal
harmonic_spurs =   2*sin(2*pi*200*t)+1*sin(2*pi*300*t)+.5*sin(2*pi*400*t); %harmonic spurs
noise = 2*randn(size(t));   %quantization noise
y = x + noise+harmonic_spurs;     % Sinusoids plus noise plus harmonic spurs

NFFT = 2^nextpow2(L); % Next power of 2 from length of y

Y = fft(y,NFFT)/(L*amplitude);  %calculate the FFT

f = Fs/2*linspace(0,1,NFFT/2);  %make the frequency vector
figure;

dBFS = 20*log10(2*abs(Y(1:NFFT/2)));  %convert FFT to dB scale
%FFT plot of signal
plot(f,dBFS)
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('dBFS')

Right about here I get stuck and don't know how to get the signal with only the noise bins.  I
would use this formula once I had only the noise bins to get the SNR:

10*log10(sum(abs(signalbins)^2)  / sum(abs(noisebins)^2))  
The SINAD would be the same except the noisebins would contain all the harmonic spurs but just
take out the input frequency spur term.

>From these values I could calculate the THD and I have already calculated the SFDR.  

Any help or tips would be appreciated.  Thanks.