Wersja kodu z komentarzami:

clear all
close all
clc

Fs = 1e3;                                 % Częstotliwość próbkowania
T = 60;                                   % Całkowity czas
t = 0:1/Fs:T-1/Fs;                        % Wektor czasu
x = chirp(t,0,t(end),500);                % Sygnał chirp
N = length(x);                            % Ilość próbek

L_seg = 10;                               % Ilość dużych odcinków
S_seg = 4;                                % Ilość małych odcinków
TotalSeg = L_seg*S_seg;                   % Całkowita ilość odcinków

segLength=N/TotalSeg;                     % Długość pojedynczego odcinka

for i=0:(TotalSeg-1)
    xSeg((i+1),:) = x((i*segLength+1):((i+1)*segLength));
    XSeg((i+1),:) = fft(xSeg((i+1),:));   % FFT z odcinków
    Xamp((i+1),:) = abs(XSeg((i+1),:));   % Amplituda
end

% Uśrednianie:
for j=1:L_seg
    Xavg(j,:) = sum(Xamp((j-1)*S_seg+1:j*S_seg,:))/S_seg;
end

df = Fs/segLength;                        % Rozdzielczość spektralna
f_vec = (0:segLength-1)*df;               % Wektor częstotliwości
t_vec = 0:6:59;
[F, TT] = meshgrid(f_vec, t_vec);

figure(1);
waterfall(F, TT, Xavg);
xlabel('Częstotliwość [Hz]');
ylabel('Czas');
zlabel('Amplituda');
xlim([0, Fs/2])

Wersja kodu bez komentarzy:

clear all
close all
clc

Fs = 1e3;
T = 60;
dt = 1/Fs;
t = 0:dt:T-dt;
x = chirp(t,0,t(end),500);
N = length(x);

L_seg = 10;
S_seg = 4;
TotalSeg = L_seg*S_seg;

segLength=N/TotalSeg;

for i=0:(TotalSeg-1)
    xSeg((i+1),:) = x((i*segLength+1):((i+1)*segLength));
    XSeg((i+1),:) = fft(xSeg((i+1),:));
    Xamp((i+1),:) = abs(XSeg((i+1),:));
end

for j=1:L_seg
    Xavg(j,:) = sum(Xamp((j-1)*S_seg+1:j*S_seg,:))/S_seg;
end

df = Fs/segLength;
f_vec = (0:segLength-1)*df;
t_vec = 0:6:59;
[F, TT] = meshgrid(f_vec, t_vec);

figure(1);
waterfall(F, TT, Xavg);
xlabel('Częstotliwość [Hz]');
ylabel('Czas');
zlabel('Amplituda');
xlim([0, Fs/2])