How to do 1D signal analysis?

I just need an example:
%% =========================
% [2.1] FFT calculation & plotting
% [2.5] Band power computation
% [1.1] Numerical integration (trapz)
% Custom PSD function (FFT-based periodogram)
% =========================
function [f, PSD] = myPSD(x, Fs)
x = x(:); % force column vector
x = x – mean(x); % remove DC offset

N = length(x);
X = fft(x);

K = floor(N/2); % one-sided spectrum length index
X = X(1:K+1);

% One-sided PSD estimate
PSD = (1/(Fs*N)) * abs(X).^2;

% Double interior bins (not DC / Nyquist)
if K > 1
PSD(2:end-1) = 2*PSD(2:end-1);
end

% Frequency vector (matches PSD length exactly)
f = (0:K)’ * (Fs/N);
end

%% Compute PSDs for original signals
[f1, PSD1] = myPSD(sig1, Fs);
[f2, PSD2] = myPSD(sig2, Fs);

%% Plot PSDs (same axes is what the lab asks)
% NOTE: If plotting linear PSD, ylabel should NOT be dB/Hz.
figure;
plot(f1, PSD1, ‘LineWidth’, 1.1); hold on;
plot(f2, PSD2, ‘LineWidth’, 1.1);
grid on;
xlim([0 200]);
xlabel("Frequency (Hz)");
ylabel("Power/Frequency (linear)");
title("PSD of Signal 1 and Signal 2");
legend("sig1","sig2");

% Optional dB version (more readable visually)
figure;
plot(f1, pow2db(PSD1), ‘LineWidth’, 1.1); hold on;
plot(f2, pow2db(PSD2), ‘LineWidth’, 1.1);
grid on;
xlim([0 200]);
xlabel("Frequency (Hz)");
ylabel("Power/Frequency (dB/Hz)");
title("PSD of Signal 1 and Signal 2 (dB)");

%% 20-100 Hz Bandpower using trapz
% [2.5] Band power computation
% [1.1] Numerical integration (trapz)
idx1 = (f1 >= 20) & (f1 <= 100);
idx2 = (f2 >= 20) & (f2 <= 100);

bandpower1 = trapz(f1(idx1), PSD1(idx1));
bandpower2 = trapz(f2(idx2), PSD2(idx2));

fprintf(‘Original bandpower (20-100 Hz): sig1 = %.6g, sig2 = %.6gn’, bandpower1, bandpower2);

%% =========================
% [1.5] Down-sampling
% Anti-aliasing before downsampling (Butterworth lowpass)
% =========================
Fs_down = 200; % new sampling rate (Nyquist = 100 Hz)
downFactor = Fs / Fs_down; % should be integer for downsample()

% Anti-alias lowpass filter (remove content above new Nyquist)
fc = 90; % use < 100 Hz for safety
order = 4;
Wn = fc / (Fs/2); % normalized cutoff

[b, a] = butter(order, Wn, ‘low’);

sig1_filt = filtfilt(b, a, sig1); % zero-phase filtering
sig2_filt = filtfilt(b, a, sig2);

% Downsample filtered signals
sig1_down = downsample(sig1_filt, downFactor);
sig2_down = downsample(sig2_filt, downFactor);

%% PSDs of anti-aliased + downsampled signals
[f1_down, PSD1_down] = myPSD(sig1_down, Fs_down);
[f2_down, PSD2_down] = myPSD(sig2_down, Fs_down);

figure;
plot(f1_down, PSD1_down, ‘LineWidth’, 1.1); hold on;
plot(f2_down, PSD2_down, ‘LineWidth’, 1.1);
grid on;
xlim([0 100]);
xlabel("Frequency (Hz)");
ylabel("Power/Frequency (linear)");
title("PSD of Downsampled Signals");
legend("sig1 down","sig2 down");

%% Bandpower (20-100 Hz) after anti-aliasing + downsampling
idx1_down = (f1_down >= 20) & (f1_down <= 100);
idx2_down = (f2_down >= 20) & (f2_down <= 100);

bandpower1_down = trapz(f1_down(idx1_down), PSD1_down(idx1_down));
bandpower2_down = trapz(f2_down(idx2_down), PSD2_down(idx2_down));

fprintf(‘Downsampled+AA bandpower (20-100 Hz): sig1 = %.6g, sig2 = %.6gn’, bandpower1_down, bandpower2_down);

%% =========================
% [1.4] Signal rectification & enveloping
% =========================
% Lab asks: rectify, then lowpass at 2 Hz to create heartbeat envelope
sig1_rect = abs(sig1_down);
sig2_rect = abs(sig2_down);

fc_env = 2; % envelope lowpass cutoff
order_env = 4;
Wn_env = fc_env / (Fs_down/2);

[b_env, a_env] = butter(order_env, Wn_env, ‘low’);

env1 = filtfilt(b_env, a_env, sig1_rect);
env2 = filtfilt(b_env, a_env, sig2_rect);

t1 = (0:length(env1)-1)’/Fs_down;
t2 = (0:length(env2)-1)’/Fs_down;

figure;
plot(t1, env1); grid on;
xlabel("Time (s)"); ylabel("Envelope amplitude");
title("Envelope of Signal 1");

figure;
plot(t2, env2); grid on;
xlabel("Time (s)"); ylabel("Envelope amplitude");
title("Envelope of Signal 2");

%% Rescale envelopes to [0,1] (helps peak detection)
env1_scaled = rescale(env1);
env2_scaled = rescale(env2);

figure;
plot(t1, env1_scaled); grid on;
xlabel("Time (s)"); ylabel("Scaled envelope");
title("Scaled Envelope of Signal 1");

figure;
plot(t2, env2_scaled); grid on;
xlabel("Time (s)"); ylabel("Scaled envelope");
title("Scaled Envelope of Signal 2");

%% =========================
% [1.2] Peak identification
% =========================
[pks1, locs1] = findpeaks(env1_scaled, ‘MinPeakProminence’, 0.1);
[pks2, locs2] = findpeaks(env2_scaled, ‘MinPeakProminence’, 0.1);

% Optional visualization of peaks
figure;
plot(t1, env1_scaled); hold on; grid on;
plot(locs1/Fs_down, pks1, ‘o’);
xlabel("Time (s)"); ylabel("Scaled envelope");
title("Signal 1 Peaks");

figure;
plot(t2, env2_scaled); hold on; grid on;
plot(locs2/Fs_down, pks2, ‘o’);
xlabel("Time (s)"); ylabel("Scaled envelope");
title("Signal 2 Peaks");

%% =========================
% [1.2] Peak timing -> Heart rate + beat-to-beat timing
% =========================
tpeaks1 = locs1 / Fs_down; % seconds
tpeaks2 = locs2 / Fs_down;

IBI1 = diff(tpeaks1); % inter-beat intervals (s)
IBI2 = diff(tpeaks2);

HR1_bpm = 60 / mean(IBI1); % average heart rate
HR2_bpm = 60 / mean(IBI2);

% Lab also asks for beat-to-beat timing variance
IBIvar1 = var(IBI1);
IBIvar2 = var(IBI2);

fprintf(‘HR sig1 = %.2f bpm, HR sig2 = %.2f bpmn’, HR1_bpm, HR2_bpm);
fprintf(‘Beat-to-beat variance: sig1 = %.6g s^2, sig2 = %.6g s^2n’, IBIvar1, IBIvar2);

%% =========================
% [1.5] Interpolation (recreate original sig1 from downsampled)
% Compare methods using RMSE
% =========================
t_down = (0:length(sig1_down)-1)’ / Fs_down;
t_orig = (0:length(sig1)-1)’ / Fs;

% Try a few methods (lab says experiment with interp1 methods)
methods = ["nearest","linear","spline","pchip","cubic"];
rmseVals = zeros(size(methods));

for k = 1:length(methods)
sig1_int = interp1(t_down, sig1_down, t_orig, methods(k))’;
rmseVals(k) = sqrt(mean((sig1 – sig1_int).^2, ‘omitnan’));
end

% Show RMSEs
disp(table(methods’, rmseVals’, ‘VariableNames’, {‘Method’,’RMSE’}));

% Pick best method (lowest RMSE)
[bestRMSE, idxBest] = min(rmseVals);
bestMethod = methods(idxBest);

% Recompute best interpolation for plotting
sig1_best = interp1(t_down, sig1_down, t_orig, bestMethod)’;

fprintf(‘Best interpolation method: %s (RMSE = %.6g)n’, bestMethod, bestRMSE);

% Plot original vs best interpolated (short window for visibility)
figure;
plot(t_orig, sig1, ‘LineWidth’, 1); hold on; grid on;
plot(t_orig, sig1_best, ‘–‘, ‘LineWidth’, 1);
xlim([0 2]); % zoom into first 2 s so differences are visible
xlabel("Time (s)"); ylabel("Amplitude");
title("Original sig1 vs Interpolated sig1");
legend("Original","Interpolated");I just need an example:
%% =========================
% [2.1] FFT calculation & plotting
% [2.5] Band power computation
% [1.1] Numerical integration (trapz)
% Custom PSD function (FFT-based periodogram)
% =========================
function [f, PSD] = myPSD(x, Fs)
x = x(:); % force column vector
x = x – mean(x); % remove DC offset