How to plot PSD in dBFS?

Question

The following code I have written measures PSD of the signal y by taking FFT first.

Fsig=26e3;           % signal frequency
Fs = 1e6;            % Sampling frequency                    
T = 1/Fs;             % Sampling period       
N = 2^15;             % Length of signal
t = (0:N-1)*T;        % Time vector

y=sin(2*pi*Fsig*t) + harmonics; // a sine function with harmonics

% find FFT+PSD
xdft = fft(y);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/N:Fs/2;
psdx_log=10*log10(psdx);

% plot PSD
figure(1)
plot(frq,psdx_log,'r')
ylim([-150 0])
grid on
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

Here's the output:

As you see the spectrum contains various harmonics in addition to the fundamental frequency. However, this plot is not ideally what I want. I want to plot PSD in dBFS where the fundamental tone, the sine wave, is at 0 dBFS so that all harmonics are measured relative to this 0 dBFS reference. What to do?