I am trying to implement envelope tracking and pulse width modulation.
In this example, I have a 1 MHz sine wave. I sample this sine wave at 3MHz and each sample will determine a duty cycle for a pulse width signal with period of 0.33us.

Pulse

sine wave

I am expecting to see the 1MHz frequency when I take the FFT of Pulse signal. but it looks like this
freq domain

what am i doing wrong?

here is my code:

pkg load signal

clear all;
close all;
% Parameters
f_sine = 1e6;        % Frequency of the sine wave (1 MHz)
t_sine = 2e-6;       % Time duration to plot (2 us)
f_sample = 3e6;      % Sampling frequency (3 MHz)
t_sample = 0:1/f_sample:t_sine;  % Time vector for sampled signal
t_cont = linspace(0, t_sine, 1000);

% Sine wave
sine_wave_cont = 1 + sin(2 * pi * f_sine * linspace(0, t_sine, 1000));  % Continuous sine wave from 0 to 2
sine_wave_sample = 1 + sin(2 * pi * f_sine * t_sample);  % Sampled sine wave

% Compute percentage of sampled amplitude compared to max amplitude
max_amplitude = max(sine_wave_sample);
percentage_amplitude = (sine_wave_sample / max_amplitude) * 100;

% Square wave parameters
T_square = 0.33e-6;  % Period of the square wave (0.33 us)

% Initialize the array for the square wave with varying duty cycles
square_wave_varying = [];

% Generate the square wave with varying duty cycle
for i = 1:length(percentage_amplitude)
    duty_cycle = percentage_amplitude(i)*0.9;
    t_current = linspace(0, T_square, 1000);   % Time vector for the current period
    if length(t_current) > 1
        square_wave_current = 0.5 * (square(2 * pi * (1/T_square) * t_current, duty_cycle) + 1);
        square_wave_varying = [square_wave_varying, square_wave_current(1:end-1)];  % Avoid overlapping periods
    end
end

% Adjust time vector for the concatenated square wave
t_varying = linspace(0, length(square_wave_varying)/f_sample, length(square_wave_varying));

% Plot continuous sine wave
figure;
plot(linspace(0, t_sine, 1000) * 1e6, sine_wave_cont);  % Time in us for x-axis
ylim([-0.5 2.5]);
xlabel('Time (us)');
ylabel('Amplitude');
title('Sine Wave at 1 MHz');
grid on;  % Turn on grid

% Plot square wave with varying duty cycle
figure;
plot(t_varying * 1e3, square_wave_varying);  % Time in us for x-axis
ylim([-0.2 1.4]);
xlabel('Time (us)');
ylabel('Amplitude');
title('Square Wave with Varying Duty Cycle');
grid on;  % Turn on grid

% Plot percentage amplitude
figure;
stem(t_sample * 1e6, percentage_amplitude, 'filled');  % Time in us for x-axis
ylim([0 110]);
xlabel('Time (us)');
ylabel('Percentage Amplitude');
title('Percentage Amplitude of Sampled Sine Wave at 3 MHz');
grid on;  % Turn on grid

% Frequency domain analysis
N = length(square_wave_varying);  % Number of points
f = (0:N-1)*(f_sample/N);         % Frequency range
Y = fft(square_wave_varying);     % Compute FFT

% Plot frequency spectrum
figure;
plot(f(1:N/2)/1e6, abs(Y(1:N/2)));  % Plot single-sided amplitude spectrum
xlabel('Frequency (MHz)');
ylabel('Magnitude');
title('Frequency Spectrum of Square Wave with Varying Duty Cycle');
grid on;  % Turn on grid