npt=8184;
ix=[0:npt-1]';
nf=396;
subh=11;  % or 33
npt2=npt/subh;
win=hanning(npt2);
win=win/sum(win);  % not necessarily right
setp_cm = zeros(npt2,nf);  setp_df = setp_cm;
seta_cm = zeros(npt2,nf);  seta_df = seta_cm;
lo=exp(j*ix*2*pi*7/33);
for fx = [1:nf]
  a=load(sprintf('raw_%3.3d.dat',fx));
  y1=a(1:npt,1);
  y2=a(1:npt,2);
  y1r=y1.*lo;   y1r=mean(reshape(y1r,subh,npt/subh)',2);
  y2r=y2.*lo;   y2r=mean(reshape(y2r,subh,npt/subh)',2);
  % common-mode
  pha_cm = (unwrap(angle(y1r))+unwrap(angle(y2r)))/2;
  pha_cm = pha_cm - mean(pha_cm);
  amp_cm = (abs(y1r)+abs(y2r))/2;
  amp_cm = log(amp_cm/mean(amp_cm));
  setp_cm(:,fx) = fft(pha_cm.*win);
  seta_cm(:,fx) = fft(amp_cm.*win);
  % differential
  pha_df = (unwrap(angle(y1r))-unwrap(angle(y2r)));
  pha_df = pha_df - mean(pha_df);
  amp_df = abs(y1r)./abs(y2r);
  amp_df = log(amp_df/mean(amp_df));
  setp_df(:,fx) = fft(pha_df.*win);
  seta_df(:,fx) = fft(amp_df.*win);
end
clk=1320e6/14;
f=clk*[0:npt2-1]'/npt;
fbin = clk/npt;
avg_cm_pha = 10*log10(mean(abs(setp_cm).^2,2)/fbin);
avg_cm_amp = 10*log10(mean(abs(seta_cm).^2,2)/fbin);
avg_df_pha = 10*log10(mean(abs(setp_df).^2,2)/fbin);
avg_df_amp = 10*log10(mean(abs(seta_df).^2,2)/fbin);
semilogx(f,avg_cm_pha, f,avg_cm_amp, f,avg_df_pha, f,avg_df_amp)
legend('common-mode phase','common-mode amplitude','differential phase','differential amplitude')
legend('location','northwest')
title('Noise Power Density')
xlim(clk*[3/npt 1/subh/2])
xlabel('f (Hz)')
ylabel('dB rad^2/Hz')
ylim([-158 -112])

% The spur at 2.8571 MHz is an artifact of the signal processing,
% the result of mixing the DC ADC offset with the digital LO.
% Other spurs at n * 0.71429 MHz are real (imperfections in DAC),
% but not interesting.

% Defining a radian to be imag(ln(complex number)) can generalize
% to amplitude units by using real(ln(complex number)).