% nco_sweep.m
%
% 4/17/2011 D. W. Hawkins (dwh@ovro.caltech.edu)
%
% NCO output frequency sweep simulation.
%
% This m-file sweeps an NCO output frequency across its
% full range and displays the output power spectrum at
% each output frequency. The spurious free dynamic range
% (SFDR) is displayed at each frequency (green line) and
% the worst SFDR is plotted in red. Each time the current
% worst-case SFDR is exceeded, the plot is held until you
% hit enter. This allows you to see what the spectrum
% looks like.
%
% The SFDR is around 5.5dB/bit of ROM address, i.e., higher
% SFDR requires more ROM bits (more angular precision).
%
% If you want 8-bits of dynamic range, ~6.02*8 = 48dB,
% out of a demodulation, then the NCO needs to have
% 48/5 ~ 9-bits of angular precision. Running this sweep
% test with [16, 10, 10] gave a SFDR of 56.24dB (5.6dB/bit).
% Even with an output of 8-bits, i.e., [16, 10, 8], the
% SFDR was still 54.97dB, so not too bad for requantization
% to 8-bits after demodulation. For [16, 8, 10], SFDR was
% 44.2dB (5.5dB/bit). For [16,6,10] SFDR was 32.3dB (5.4dB/bit).
%
% Note that Altera's NCO compiler v10.1 will no longer let
% you select lower than 10-bits of output magnitude precision
% (at least on the Stratix IV GX test I was doing).
%

% Bit-widths
B_acc  = 16;
B_addr = 8; 
B_data = 10;

B_acc  = 8;
B_addr = 8; 
B_data = 8;

% Number of samples per spectra
Nt = 2048;

% NCO clock frequency
f_clk = 100e6;

% Window function
w = kaiser(Nt,4*pi)';
wcg = mean(w);

% Test frequencies
f_out = linspace(1,49,Nt)*1e6;
%f_out = [1:0.25:49]*1e6;
%f_out = linspace(1,49,1234)*1e6;

% NCO frequency step size (based on accumulator size)
f_step = f_clk/2^B_acc;

% Sample index
n = [0:Nt-1]; 

% Frequency channels (FFT output format)
f = [0:Nt-1]/Nt*f_clk; 

% Cos and Sin ROMs
% * 2^B_addr locations quantized to B_data bits
rom_addr = [0:2^B_addr-1];
cos_rom = round(cos(2*pi*rom_addr/2^B_addr)*2^(B_data-1));
sin_rom = round(sin(2*pi*rom_addr/2^B_addr)*2^(B_data-1));

% Generate a figure at each frequency
M = length(f_out);
sfdr_dB = -500;
stop_plot = 0;
for m = 1:M,

	% Convert requested f_out to an FTW
	ftw = round(f_out(m)/f_step);
	
	% Actual output frequency
	f_act = ftw*f_step;
	
	fprintf('%4d: f_out = %6.3f MHz, f_act = %6.3f MHz\n', m, f_out(m)/1e6, f_act/1e6)
	
	% Accumulator output at each time instant
	acc_out = ftw*n;
	
	% Accumulator-to-ROM addressing
	%  * Truncate the accumulator output to ROM address units
	%  * Wrap-around the ROM addresses
	rom_addr = mod(floor(acc_out/2^B_acc*2^B_addr),2^B_addr);

	% 'ROM' cos and sin lookup
	% * add 1, due to MATLAB's indexing scheme
	x_cos = cos_rom(rom_addr+1);
	x_sin = sin_rom(rom_addr+1);

	% Spectra
	% * normalize the signal by 2^(B_data-1) to put it in
	%   fractional integer format (-1.0 to 1.0)
	% * normalize by (Nt/2) for the expected magnitude of
	%   the cosine/sine tranform (0dB peak)
	% * normalize by the window coherent gain
	%
	X_cos = fft(w.*x_cos)/2^(B_data-1)*2/Nt/wcg;
	X_sin = fft(w.*x_sin)/2^(B_data-1)*2/Nt/wcg;
	
	% Power spectra
	Rxx_cos = abs(X_cos).^2;
	Rxx_sin = abs(X_sin).^2;
	
	% Find the worst-case harmonic peak
	% * mask off the channels near the signal
	N_act = round(f_act/f_clk*Nt);
	mask = ones(1,Nt);
	mask(     N_act + [-5:5]) = 0;
	mask(Nt - N_act + [-5:5]) = 0;
	mask([[1:3] [Nt-2:Nt]])  = 0; % DC
	mask(Nt/2 + [-3:3]) = 0; % Nyquist
	max_cos = max(Rxx_cos.*mask);
	max_sin = max(Rxx_sin.*mask);
	max_dB  = 10*log10(max([max_cos max_sin]));
	
	if (max_dB > sfdr_dB)
		sfdr_dB = max_dB;
		if (m > 1)
			stop_plot = 1;
		end
	end

	figure(1)
	hold off
	plot(f(1:Nt/2+1)/1e6,10*log10(Rxx_cos(1:Nt/2+1)))
	hold on
	plot(f(1:Nt/2+1)/1e6,10*log10(Rxx_sin(1:Nt/2+1)),'r')
	plot([f(1) f(Nt/2+1)]/1e6,[0 0],'k:')
	plot([f(1) f(Nt/2+1)]/1e6,[1 1]*max_dB,'g:')
	plot([f(1) f(Nt/2+1)]/1e6,[1 1]*sfdr_dB,'r:')
	axis([[f(1) f(Nt/2+1)]/1e6 -100 5])
	xlabel('Frequency (MHz)')
	ylabel('Power (dB)')
	drawnow

	if (stop_plot == 1)
		% Stop each time a new maximum occurs to
		% see where those frequencies are
		% * after about 10MHz the max does not get
		%   exceeded again (or too often)
		input('Hit enter','s');
		stop_plot = 0;
	end	
end
fprintf('Worst-case SFDR is %.2f dB\n', -sfdr_dB)
fprintf('SFDR/ROM width %.2f dB per angle bit\n', -sfdr_dB/B_addr)