In the large scale, the modulated spectrum looks actually plaussible up to ~500MHz.
Don't know what's the "real" sample rate? 1200 or 2400 MSa/s?
If it is 1200, then the cut-off beyond 500 MHz is to be expected, and about -34dBc at 500 Mhz is expected too, for a 10 MHz square wave.
The unmodulated spectrum does not look like a square wave spectrum, but it lacks high-frequency contents.
Were different rise/fall times configured for the modulated and unmodulated test case?
Maybe you can check the waveforms and the rise times with a scope?
There are even carrier harmonics present, which are unexpected for a square wave. However their amplitudes are relatively low.
Please find attached the spectrum I would expect for a 10 MHz square wave, properly band-limited to 480 MHz (=> -120dBc at 600 MHz), and then PM-modulated with 1° deviation and a 1kHz sine wave signal. Sorry can't use 100 Hz, as this would require 10x more samples, and Gnuplots runs out of memory then. For the same reason I also can't use a smaller RBW, therefore the lobes in the zoom-in images are pretty wide. But each lobe represents just a single frequency line, not a spread spectrum.
EDIT: Updated images. Noticed that the given frequency plan allows a rectangular window w/o introducing leakage. This enables a narrower RBW for the spectrum plot.
EDIT: Added script which generates the plot
EDIT: 2023-04-27 Did a couple of script enhancements
- square wave edges can be shaped according to given rise time and shaping method
- alternative band-limiting methods (FIR requires too many taps)
- minimize required FFT resolution and size
- plot the (rise/fall-time adjusted and bandwidth limited) wavetable which is then re-sampled via DDS
- find and print peaks
pkg load signal
more off
k = 1000
M = 1E6
fc = 10*M % carrier freq [Hz]
waveform = "square" % square/sine
edges = "gaussian" % gaussian/sinc/sinc2/sinc3 (only for square wave)
risetime = 0 % risetime [s] (only for square wave)
%risetime = 2e-9 % SDG6000X ?
%risetime = 8.4e-9 % SDG2000X ?
modulation = "PM" % FM/PM/none
fmod = 1*k % modulating signal (sine) frequency [Hz]
df = 1 % FM frequency deviation [Hz]
dp = 1 % PM phase deviation [degrees]
fs = 1200*M % DDS sample rate [Sa/s]
%fs = 300*M % DDS sample rate [Sa/s]
bandlimit = "RC" % RC/brickwall/FIR
BW = 0.4*fs % must be < 0.5*fs (for RC and FIR)
ndds = 4*1024*1024 % wavetable size
%-------------------------------------------------------------
% Find reasonable binsize and fftsize
% Try to use wider bins and rectangular window if possible
% without leakage, use narrower bins and Kaiser window otherwise.
% Quantiize frequencies to fix floating point inaccuracies.
fs = floor(fs*1024+0.5)/1024;
fc = floor(fc*1024+0.5)/1024;
fmod = floor(fmod*1024+0.5)/1024;
if mod(2*fc,fmod) == 0 && mod(2*fs,fmod) == 0
binsize = fmod/2
fftsize = 2*fs/fmod
window = "rect"
RBW_3dB = 0.88*binsize
else
fftsize = floor(10*fs/fmod+0.5)
binsize = fs/fftsize
window = "kaiser"
RBW_3dB = 2.24*binsize
end
% waveform prototype, amplitude 1 units peak
if strcmp(waveform,"square")
wavetable = [-ones(1,ndds/2) ones(1,ndds/2)] * sqrt(0.5);
elseif strcmp(waveform,"sine")
wavetable = sin([0:ndds-1]*(2*pi/ndds));
else
error ("Illegal waveform");
end
% edge shaping for square wave
shaper = [ 1 zeros(1,ndds-1) ];
if strcmp(waveform,"square") && risetime > 0
if strcmp(edges,"gaussian")
% Gaussian shaper (80% center quantile is 1.2816*sigma)
shaper = ifftshift(gausswin(ndds,1.2816/fc/risetime)');
elseif strcmp(edges,"sinc")
% sinc shaper (moving average)
shaper = ones(1,floor(1.25*risetime*fc*ndds));
shaper = [ shaper zeros(1,ndds-length(shaper)) ];
elseif strcmp(edges,"sinc2")
% sinc^2 shaper
shaper = ones(1,floor(0.905*risetime*fc*ndds));
shaper = fftconv(shaper,shaper);
shaper = [ shaper zeros(1,ndds-length(shaper)) ];
elseif strcmp(edges,"sinc3")
% sinc^3 shaper
shaper = ones(1,floor(0.76*risetime*fc*ndds));
shaper = fftconv(fftconv(shaper,shaper),shaper);
shaper = [ shaper zeros(1,ndds-length(shaper)) ];
else
error ("Illegal shaper");
end
shaper /= sum(shaper);
else
shaper = [ 1 zeros(1,ndds-1) ];
end
% band-limit samples in wavetable and apply shaper
if strcmp(bandlimit,"RC")
% apply raised cosine roll-off BW...0.5*fs
low = BW/fc;
high = 0.5*fs/fc;
idx = [ceil(low):1:floor(high)];
WT = fft(wavetable).*fft(shaper);
WT(idx) .*= cos((idx-low)/(high-low)*pi/2);
WT(ceil(high)+1:end-ceil(high)) = 0;
WT(ndds/2+2:end) = conj(fliplr(WT(2:ndds/2)));
wavetable = real(ifft(WT));
elseif strcmp(bandlimit,"brickwall")
% apply ideal brickwall at Nyquist
keep = floor(0.5*fs/fc);
WT = fft(wavetable).*fft(shaper);
WT(keep+1:end-keep) = 0;
wavetable = real(ifft(WT));
elseif strcmp(bandlimit,"FIR")
% FIR filter with BW...0.5*fs transition
% band and 120dB stopband attenuation
[n,Wn,beta,ftype] = kaiserord([BW/ndds/fc 0.5*fs/ndds/fc], [1 0], [0.001 0.000001], 1);
if n+1>ndds
error("Too many filter taps required, choose a different 'bandlimit' method")
end
taps = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
taps = [ taps zeros(1,ndds-length(taps)) ];
wavetable = real(ifft(fft(wavetable).*fft(taps).*fft(shaper)));
else
error("Illegal bandlimiting method")
end
% modulation signal (sine)
mod_signal = sin(2*pi*fmod*[0:fs/fmod-1]/fs);
mod_signal = repmat(mod_signal,1,ceil(fftsize/length(mod_signal)))(1:fftsize);
% generate modulated carrier via DDS
if strcmp(modulation,"FM")
phases = cumsum(ndds/fs*(fc+mod_signal*df));
elseif strcmp(modulation,"PM")
phases = [0:length(mod_signal)-1]*ndds/fs*fc + mod_signal*ndds*dp/360;
elseif strcmp(modulation,"none")
phases = [0:length(mod_signal)-1]*(fs*fc/ndds);
else
error("Illegal modulation")
end
signal = wavetable(1+floor(mod(phases + ndds, ndds)));
% 16-bit quantization
%signal = floor(signal * 32768 + 0.5) / 32768;
% 14-bit quantization
%signal = floor(signal * 8192 + 0.5) / 8192;
% 14-bit quantization + random noise
%signal = floor(signal * 8192 + 0.5 + 0.5*randn(1,length(signal))) / 8192;
% Window function for FFT
if strcmp(window,"kaiser")
w = kaiser(length(signal)+1,17)(1:end-1)';
else
w = ones(1, length(signal));
end
w /= mean(w);
% plot spectrum
frequencies = [0:length(signal)-1]*fs/length(signal);
magnitude = abs(fft(signal.*w)/length(signal)*2);
figure 1
clf
plot(frequencies/M,20*log10(magnitude+1e-12))
grid on
title(sprintf("%s, deviation %g°, fc=%gMHz (square), fmod=%gkHz (sine)", modulation, dp, fc/M, fmod/k))
xlabel("MHz")
ylim([-100 5])
xlim([0 0.5*fs/M])
%xlim([10-0.007 10+0.007]) % zoom-in to 10MHz
%xlim([450-0.007 450+0.007]) % zoom-in 450 MHz
% plot wavetable
figure 2
clf
plot([0:ndds-1]/(fc*ndds*1e-9),wavetable)
title("Wavetable")
xlabel("ns")
grid on
% find and print peaks > -120dB
idx = 2:length(signal)/2-1;
peaks = find(frequencies(idx)<=BW &
magnitude(idx)>1e-6 &
magnitude(idx)>magnitude(idx-1) &
magnitude(idx)>magnitude(idx+1)) + 1;
for f = peaks
printf("peak: %10.6f MHz: %6.1f dB\n",frequencies(f)/M,20*log10(magnitude(f)+1e-12));
end
Home
Help
Search
Login
Register
