Can anyone help me with my code? I want to apply the windowing function to my code with a duration of 200 ms to observe the leakage.

function [y, signalT, notch_intensity] = signalfunccycle(Gridlength, startfreq, endfreq, notchstart1, notchend1, datalength)
close all;
clear all;
PIx2 = 2 * pi; % used to calculate the frequency values of the output signal
fstart = 0;
Gridlength = 200e-3;
startfreq = 40e3;
endfreq = 200e3;
stimeStep = 5e-8;
datalength = round(Gridlength / stimeStep);

notchstart1 = 116000;
notchend1 = 119000;
dataN = datalength;
Stimelength = Gridlength;
fend = 10^6;
Fs = dataN / Gridlength;
NFFT = dataN; % Number of points used in the FFT analysis
fstep = (Fs – fstart) / NFFT; % Point step

T = 1 / Fs;
Stime = linspace(0, Stimelength, dataN);
t = (0:dataN-1) * T; % Length of time
signal = zeros(1, length(t));
f1 = notchstart1;
f2 = notchend1;
f3 = notchstart1;
f4 = notchend1;

for freq = startfreq:0.5*10^3:endfreq
if (f1 <= freq && freq <= f2) || (f3 <= freq && freq <= f4)
continue;
end
signal = sin(PIx2 * freq * t + pi) + signal;
end

%% Time-domain signal graphs
figure(1);
plot(Stime, signal);
grid on;
title(‘signal’);
xlabel(‘time’);
ylabel(‘amplitude’);

%% FFT
fMax = NFFT / 2 + 1; % The maximum frequency to plot
signalFFT0 = fft(signal, NFFT); % FFT of the signal
signalFFT_phase = angle(fft(signal, NFFT)); % Phase of the FFT result
signalFFT = abs(fft(signal, NFFT)); % Absolute value (modulus) of the FFT result
signalFFTShow = 2 * signalFFT / dataN; % Actual amplitude values
signalFFTShow(1) = signalFFTShow(1) / 2; % DC component
f = Fs / 2 * linspace(0, 1, fMax); % Actual frequencies

% Plot FFT of the signal
figure(2);
plot(f, signalFFTShow(1:fMax), ‘r’, ‘LineWidth’, 3); % Red color, line width set to 3
grid on; % Display grid lines
title(‘FFT Signal’);
xlabel(‘Frequency’);
ylabel(‘Amplitude’);

%% Phase function modulation
signalstart = startfreq;
signalend = endfreq;
Istart = round(signalstart / fstep);
Iend = round(signalend / fstep);
Itotal = Iend – Istart;

N1start1 = round(f1);
N1end1 = round(f2);
N1start2 = round(f3);
N1end2 = round(f4);
signalphase(1:NFFT) = 0;

for I = 1:NFFT
if (I <= N1end1 && I >= N1start1) || (I <= N1end2 && I >= N1start2)
signalphase(I) = 0;
else
signalphase(I) = pi / 2 * (I – Istart) + 1 * pi / Itotal / 2 * (I – Istart)^2;
end
end

signalF = 1/2 * signalFFTShow .* exp(1i * (signalphase + signalFFT_phase));
signalT = ifft(NFFT * signalF);
t1 = (0:NFFT-1) * T;

%% Phase spectrogram
x1 = 1:NFFT;
figure(3);
plot(x1, signalphase)
xlabel(‘frequency’);
ylabel(‘Phase’);
title(‘Phase spectrum’);

%% Signal waveform after modulation
signaltime = (0:NFFT-1) * T;
figure(4);
plot(signaltime, real(signalT));
xlabel(‘time’);
ylabel(‘amplitude’);
signalT_real = real(signalT);
savefile = ‘signalT_real.txt’;
x = Stime’;
y = signalT_real’;
save(savefile, ‘x’, ‘y’, ‘-ASCII’)

%% STFT Calculation
NN = 10;
win_size = round(datalength / NN);
overlap = 0;
fft_size = win_size * 2 + 1;
stft_result = [];

for i = 1:win_size*(1-overlap):length(signalT)-win_size
window = signalT(i:i+win_size-1);
fft_result = fft(window, fft_size+1);
stft_result = [stft_result; abs(fft_result(1:fft_size))];
end

time_axis = linspace(0, Gridlength, size(stft_result, 1));
freq_axis = linspace(fstart, fend, fft_size);

% Extract intensity for notch frequencies (116 to 119 kHz)
notch_indices = find(freq_axis >= 116000 & freq_axis <= 119000);
notch_intensity = mean(stft_result(:, notch_indices), 2);

% Plot the STFT spectrogram
figure(5);
imagesc(time_axis * 1000, freq_axis / 1000, stft_result.’);
axis([Gridlength * 950 Gridlength * 1050 110e3 / 1000 125e3 / 1000]);
colorbar;
ylabel(‘Frequency (kHz)’);
xlabel(‘Time (ms)’);

% % Plot the intensity over time for the notch frequencies
% figure(6);
% plot(time_axis * 1000, notch_intensity, ‘LineWidth’, 2);
% xlabel(‘Time (ms)’);
% ylabel(‘Intensity’);
% title(‘Intensity at Notch Frequencies (116 to 119 kHz)’);
% grid on;
endfunction [y, signalT, notch_intensity] = signalfunccycle(Gridlength, startfreq, endfreq, notchstart1, notchend1, datalength)
close all;
clear all;
PIx2 = 2 * pi; % used to calculate the frequency values of the output signal
fstart = 0;
Gridlength = 200e-3;
startfreq = 40e3;
endfreq = 200e3;
stimeStep = 5e-8;
datalength = round(Gridlength / stimeStep);