https://dsp.stackexchange.com/posts/76761/

simulate_qam_spectrum.m

clear('all');
close('all');
clc();

fsdac = 12e9;  % sampling rate of DAC system
fnco = 2e9;    % carrier frequency
fout = 450e3;  % I/Q baseband channel modulated frequency
tsim = 100000e-9; % simulation time

M = 16;
message = 0:M-1;
symbols = qammod(message,M,'bin', 'UnitAveragePower',true, 'PlotConstellation', true);

% oversampling factor in order to get more datapoints / simulate a zero-order-hold at the modulator output
ovs = 5;  

t = linspace(0, tsim, fsdac * tsim);
t(end) = []; % 
tovs = linspace(0, tsim, fsdac * tsim * ovs + 1 - ovs);
tovs(end) = [];

f = linspace(0, fsdac * ovs, size(tovs, 2));
iqgrid = figure();

% for each symbol
for kfo = 1 : size(symbols, 2)
    % baseband I/Q modulation
    % QAM is amplitude modulation --> take the amplitude from the symbols and apply it to the baseband frequency
    xi = real(symbols(kfo)) * cos( 2*pi*fout*t );
    xq = imag(symbols(kfo)) * -sin( 2*pi*fout*t );
    
    % NCO carrier frequency
    ncoi = cos( 2*pi*fnco*t );
    ncoq = sin( 2*pi*fnco*t );
    
    % I/Q Modulation (single channel output)
    outi = xi.*ncoi + xq.*ncoq;
    
    % Quantization of modulated signal
    % signed 16 bit mit 14 bit Nachkommastellen
    x = storedIntegerToDouble( fi( outi, true, 16, 14) ) ./ (2^14); 
    
    % digital to analog: ZOH
    xovs = repelem(x,1,ovs);
    % digital to analog: Return to Zero / Pulse
%     xovs = upsample(x, ovs);    

    % Nonlinearity (analoge amplifier after DAC)
    % tanh shows saturation effect like power amplifiers
    k = 1.0;
    yovs = tanh(xovs .* k) ./ k;

    % FFT
%     yovs(kfo, :) = 20.*log10(abs(fft(xovs) .* (2 / size(tovs, 2)))); %#ok<SAGROW>
    specx = 20.*log10(abs(fft(xovs) .* (2 / size(tovs, 2)))); 
%     specy = 20.*log10(abs(fft(yovs) .* (2 / size(tovs, 2)))); 
    
    % Coordinates on I/Q grid
    coord = 1+mod((kfo-1)*sqrt(M),M-1);
    if kfo == size(symbols,2)
        coord = size(symbols,2);
    end
    subplot(sqrt(M),sqrt(M), coord );
    grid('on');
    hold('on');
    if(coord == (M-sqrt(M)+1))
        xlabel('Frequency in Hz');
        ylabel('Power in dBc');
    end
    
    plot(f, specx);
%     plot(f, specy);

    axis([fnco - 10e6, fnco + 10e6, -150, 0]);
    titlestr = sprintf("I/Q = [%.3f, %.3f]", real(symbols(kfo)), imag(symbols(kfo)));
    title(titlestr);
end

% titlestr2 = sprintf("QAM-%d Spectrum Simulation with nonlinear Distortion", M); 

gridax=axes(iqgrid,'visible','off'); 
gridax.Title.Visible='on';
gridax.XLabel.Visible='on';
gridax.YLabel.Visible='on';
ylabel(gridax,'Quadrature Amplitude');
xlabel(gridax,'In-phase Amplitude');
titlestr2 = sprintf("QAM-%d Spectrum Simulation", M); 
sgtitle(titlestr2);

set(gcf, 'Position',  [100, 100, 800, 800])