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 = %.6g\n', 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 = %.6g\n', 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 bpm\n', HR1_bpm, HR2_bpm);
fprintf('Beat-to-beat variance: sig1 = %.6g s^2, sig2 = %.6g s^2\n', 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");

[Mathieu NOE]

what is your question ?

[Joshua]

%% =========================================================

%% 1.1 Numerical integration & differentiation

%% =========================================================

% PRACTICE QUESTION:

% Given a time vector t (s) and signal x(t), compute:

% 1) dx/dt (numerical derivative)

% 2) the area under x(t) from t = 2 s to t = 5 s

% --- CODED ANSWER ---

% Assume t and x are already given

t = t(:); % force column vector (keeps dimensions consistent)

x = x(:);

dt = mean(diff(t)); % sample interval (assuming uniform sampling)

% Why: needed to convert "change per sample" into "change per second"

% 1) Numerical derivative (same length as x)

dxdt = gradient(x, dt);

% Why: derivative tells you how fast the signal is changing (slope/rate of change)

% Important for identifying rapid transitions, motion, onset changes, etc.

% 2) Numerical integration over a time window [2, 5] s

idx = (t >= 2) & (t <= 5);

% Why: selects only the time region we care about

area_2_5 = trapz(t(idx), x(idx));

% Why: integration adds up the signal over time (area under curve)

% Important for total quantity/accumulated signal/power-like calculations

% Optional plot

figure;

plot(t, x); hold on; grid on

plot(t(idx), x(idx), 'LineWidth', 1.5)

xlabel('Time (s)'); ylabel('x(t)');

title('1.1 Integration + Differentiation');

legend('Signal','Integrated Region');

%% =========================================================

%% 1.2 Peak & trough identification

%% =========================================================

% PRACTICE QUESTION:

% Given a 1D signal x(t), identify peaks and troughs using findpeaks.

% Return peak/trough amplitudes and locations, and plot them on the signal.

% --- CODED ANSWER ---

t = t(:);

x = x(:);

% Find peaks (tune parameters if needed)

[pks, locs_pk] = findpeaks(x, t, 'MinPeakProminence', 0.1);

% Why: finds meaningful local maxima (ignores tiny noisy bumps)

% Important for heart beats, pulse peaks, EMG bursts, etc.

% Find troughs by flipping the signal

[trs_neg, locs_tr] = findpeaks(-x, t, 'MinPeakProminence', 0.1);

trs = -trs_neg;

% Why: MATLAB finds peaks, not troughs, so we flip the signal to reuse findpeaks

% Important when you need minima timing/amplitude too

% Plot

figure;

plot(t, x); hold on; grid on

plot(locs_pk, pks, 'o', 'LineWidth', 1.5)

plot(locs_tr, trs, 'x', 'LineWidth', 1.5)

xlabel('Time (s)'); ylabel('Signal');

title('1.2 Peaks and Troughs');

legend('x(t)','Peaks','Troughs');

%% =========================================================

%% 1.3 Signal quality assessment

%% =========================================================

% PRACTICE QUESTION:

% Given a raw signal x_raw and a filtered version x_filt, estimate SNR (dB)

% assuming:

% signal = x_filt

% noise = demeaned raw - filtered

% --- CODED ANSWER ---

x_raw = x_raw(:);

x_filt = x_filt(:);

% Remove DC from raw signal first

x_dm = x_raw - mean(x_raw, 'omitnan');

% Why: DC offset is a constant shift, not useful signal content

% Important because offsets can distort power/SNR and filtering steps

% Estimate noise as difference between demeaned raw and filtered signal

noise_est = x_dm - x_filt;

% Why: what's removed by filtering is treated as "noise"

% Important for quantifying how much junk was in the raw signal

% Power-based SNR (works well for signal quality assessment)

P_signal = mean(x_filt.^2, 'omitnan');

P_noise = mean(noise_est.^2, 'omitnan');

% Why: power = average squared amplitude (standard way to compare signal vs noise)

SNR_dB = 10*log10(P_signal / P_noise);

% Why: expresses signal quality on a dB scale (easy to compare)

% Higher SNR = cleaner signal

fprintf('Estimated SNR = %.2f dB\n', SNR_dB);

% Optional plot

figure;

plot(x_dm); hold on; plot(x_filt); plot(noise_est);

grid on; title('1.3 Signal Quality Assessment');

legend('Demeaned Raw','Filtered (Signal)','Estimated Noise');

%% =========================================================

%% 1.4 Signal rectification & enveloping

%% =========================================================

% PRACTICE QUESTION:

% Given a filtered EMG/audio-type signal x and sampling frequency Fs:

% 1) Full-wave rectify the signal

% 2) Compute an envelope using a lowpass Butterworth filter

% 3) Plot raw/rectified/envelope

% --- CODED ANSWER ---

x = x(:);

% 1) Full-wave rectification

x_rect = abs(x);

% Why: makes all values positive so oscillations don't cancel out

% Important because EMG/audio oscillates around zero, but we care about "magnitude/activity"

% 2) Lowpass filter for envelope (example cutoff = 3 Hz)

fc_env = 3;

order_env = 4;

Wn_env = fc_env/(Fs/2);

% Why: lowpass keeps the slow trend and removes fast wiggles

% This smooth trend is the "envelope" (overall activation pattern)

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

x_env = filtfilt(b_env, a_env, x_rect);

% Why filtfilt: zero-phase filtering (doesn't shift peaks in time)

% Important for timing-based analysis (onsets, peaks, heart/EMG events)

% Time vector (if not provided)

t = (0:length(x)-1)'/Fs;

% 3) Plot

figure;

plot(t, x, 'DisplayName','Filtered signal'); hold on; grid on

plot(t, x_rect, 'DisplayName','Rectified');

plot(t, x_env, 'LineWidth',1.5, 'DisplayName','Envelope');

xlabel('Time (s)'); ylabel('Amplitude');

title('1.4 Rectification and Envelope');

legend;

%% =========================================================

%% 1.5 Interpolation & down-sampling

%% =========================================================

% PRACTICE QUESTION:

% Downsample a signal from Fs to 200 Hz (with anti-aliasing), then interpolate

% it back to the original time vector and compute RMSE.

% --- CODED ANSWER ---

x = x(:);

Fs_new = 200;

D = Fs / Fs_new;

% Why: downsample factor tells MATLAB how many samples to skip

% Important for reducing data size and processing load

% Anti-alias filter BEFORE downsampling (lowpass below new Nyquist)

fc_aa = 90;

[b_aa, a_aa] = butter(4, fc_aa/(Fs/2), 'low');

x_aa = filtfilt(b_aa, a_aa, x);

% Why: removes frequencies that would fold into lower frequencies (aliasing)

% Important because downsampling without this can corrupt the signal

% Downsample

x_down = downsample(x_aa, D);

% Why: actually reduces the sample rate by keeping every D-th sample

% Build time vectors

t_orig = (0:length(x)-1)'/Fs;

t_down = (0:length(x_down)-1)'/Fs_new;

% Interpolate back to original time points

x_interp = interp1(t_down, x_down, t_orig, 'cubic');

% Why: estimates missing points between downsampled samples

% Important for reconstructing/upsampling and comparing methods

% Compute reconstruction error

RMSE = sqrt(mean((x - x_interp).^2, 'omitnan'));

% Why: RMSE quantifies how close the reconstruction is to the original

% Lower RMSE = better interpolation

fprintf('Interpolation RMSE = %.6f\n', RMSE);

% Plot short segment

figure;

plot(t_orig, x); hold on; grid on

plot(t_orig, x_interp, '--');

xlim([0 2]);

xlabel('Time (s)'); ylabel('Amplitude');

title('1.5 Downsample + Interpolate');

legend('Original','Interpolated');

%% =========================================================

%% 2.1 FFT calculation & plotting

%% =========================================================

% PRACTICE QUESTION:

% Compute and plot the one-sided PSD of a signal x sampled at Fs.

% --- CODED ANSWER ---

x = x(:);

x = x - mean(x);

% Why: removes DC offset so the 0 Hz component doesn't dominate the spectrum

% Important for cleaner frequency analysis

N = length(x);

X = fft(x);

% Why: FFT converts time-domain signal into frequency-domain components

% Important for seeing what frequencies are present

K = floor(N/2);

X = X(1:K+1);

% Why: for real signals, the spectrum is symmetric, so we keep only one side

% Important for one-sided PSD plotting (0 to Nyquist)

% One-sided PSD estimate

PSD = (1/(Fs*N)) * abs(X).^2;

if K > 1

PSD(2:end-1) = 2*PSD(2:end-1);

end

% Why: converts FFT magnitude to power per Hz and accounts for removed negative frequencies

% Important for bandpower and frequency-domain comparisons

% Frequency axis (must match PSD length)

f = (0:K)' * (Fs/N);

% Why: maps each PSD bin to its actual frequency in Hz

% Plot in dB/Hz

figure;

plot(f, pow2db(PSD)); grid on

xlabel('Frequency (Hz)');

ylabel('Power/Frequency (dB/Hz)');

title('2.1 One-Sided PSD via FFT');

%% =========================================================

%% 2.2 FIR filtering

%% =========================================================

% PRACTICE QUESTION:

% Design an FIR lowpass filter with cutoff 40 Hz and apply it to x.

% Plot raw and filtered signals.

% --- CODED ANSWER ---

x = x(:);

% FIR = Finite Impulse Response (no feedback)

fir_order = 100;

fc = 40;

Wn = fc/(Fs/2);

% Why: normalized cutoff is required by MATLAB filter design functions

% Important because filter design always uses frequency relative to Nyquist

% Design FIR filter coefficients

b_fir = fir1(fir_order, Wn, 'low');

% Why: creates FIR filter coefficients for a lowpass filter

% Important for removing high-frequency noise while keeping lower frequencies

% Apply FIR filter

% filter(...) is causal; filtfilt(...) removes phase delay

x_fir = filtfilt(b_fir, 1, x);

% Why filtfilt: applies filter forward and backward to avoid time shift

% Important when timing/peak location matters

% Plot

t = (0:length(x)-1)'/Fs;

figure;

plot(t, x); hold on; grid on

plot(t, x_fir, 'LineWidth',1.2)

xlabel('Time (s)'); ylabel('Amplitude');

title('2.2 FIR Lowpass Filtering');

legend('Raw','FIR filtered');

%% =========================================================

%% 2.3 Notch filtering for mains interference removal

%% =========================================================

% PRACTICE QUESTION:

% Remove 60 Hz mains interference from a signal x sampled at Fs.

% --- CODED ANSWER ---

x = x(:);

% Mains interference is usually 60 Hz (in your course context)

f_notch = 60;

BW = 4;

% Why: notch is a narrow filter that removes one specific interference frequency

% Important because power-line noise often contaminates biosignals

% Normalize for iirnotch

wo = f_notch/(Fs/2);

Q = f_notch/BW;

bw = wo/Q;

% Why: iirnotch needs normalized center frequency and bandwidth

% Q controls how narrow/sharp the notch is

% Design and apply notch filter

[b_n, a_n] = iirnotch(wo, bw);

x_notch = filtfilt(b_n, a_n, x);

% Why: removes 60 Hz noise while keeping most of the rest of the signal

% Important for cleaner EMG/EEG and better PSD/features

% Plot (time-domain)

t = (0:length(x)-1)'/Fs;

figure;

plot(t, x); hold on; grid on

plot(t, x_notch);

xlabel('Time (s)'); ylabel('Amplitude');

title('2.3 60 Hz Notch Filtering');

legend('Raw','Notch filtered');

%% =========================================================

%% 2.4 Signal identification

%% =========================================================

% PRACTICE QUESTION:

% You are given a signal x and Fs. Use time and frequency plots to identify

% whether it likely contains low-frequency drift, mains noise (60 Hz), or

% high-frequency activity. (Return a brief interpretation.)

% --- CODED ANSWER ---

x = x(:);

t = (0:length(x)-1)'/Fs;

% Time-domain plot

figure;

plot(t, x); grid on

xlabel('Time (s)'); ylabel('Amplitude');

title('2.4 Signal Identification - Time Domain');

% Why: some artifacts are easier to see in time (drift, bursts, spikes)

% PSD (reuse FFT/PSD workflow)

x_dm = x - mean(x);

N = length(x_dm);

X = fft(x_dm);

K = floor(N/2);

X = X(1:K+1);

PSD = (1/(Fs*N))*abs(X).^2;

if K > 1

PSD(2:end-1) = 2*PSD(2:end-1);

end

f = (0:K)'*(Fs/N);

figure;

plot(f, pow2db(PSD)); grid on

xlabel('Frequency (Hz)'); ylabel('Power/Frequency (dB/Hz)');

title('2.4 Signal Identification - PSD');

% Why: frequency plot helps identify what kind of signal/artifact is present

% (e.g., 60 Hz noise, low-frequency drift, band-limited activity)

% Example "identification" checks (simple automatic flags)

% 1) Check for 60 Hz mains peak

idx60 = (f >= 58 & f <= 62);

has60Hz = any(PSD(idx60) > 5*median(PSD));

% Why: a strong narrow peak around 60 Hz often means mains interference

% 2) Check for low-frequency drift (< 2 Hz power)

idxLow = (f >= 0 & f <= 2);

idxMid = (f > 2 & f <= 20);

lowPower = trapz(f(idxLow), PSD(idxLow));

midPower = trapz(f(idxMid), PSD(idxMid));

hasDrift = lowPower > midPower;

% Why: if very low frequencies dominate, the signal likely has baseline drift

% 3) Simple printed interpretation

if has60Hz

disp('Likely mains interference near 60 Hz detected.');

end

if hasDrift

disp('Likely low-frequency drift/baseline wander present.');

end

if ~has60Hz && ~hasDrift

disp('Signal appears cleaner; inspect dominant bands manually.');

end

%% =========================================================

%% 2.5 Band power computation

%% =========================================================

% PRACTICE QUESTION:

% Given a PSD (f, PSD), compute the bandpower from 20 to 100 Hz.

% --- CODED ANSWER ---

% Assume f and PSD are already computed

f_low = 20;

f_high = 100;

idx_band = (f >= f_low) & (f <= f_high);

% Why: selects only the frequencies inside the band of interest

band_power_20_100 = trapz(f(idx_band), PSD(idx_band));

% Why: bandpower = area under the PSD in that band

% Important for comparing signal energy in specific frequency ranges

fprintf('Bandpower (20-100 Hz) = %.6f\n', band_power_20_100);

% Optional visualize selected band

figure;

plot(f, PSD); hold on; grid on

plot(f(idx_band), PSD(idx_band), 'LineWidth',1.5)

xlabel('Frequency (Hz)'); ylabel('Power/Frequency (linear)');

title('2.5 Band Power (20-100 Hz)');

legend('PSD','Integrated Band');

%% =========================================================

%% 3.1 2D image processing

%% =========================================================

% PRACTICE QUESTION:

% Load a noisy grayscale image, remove salt-and-pepper noise using medfilt2,

% compare to averaging filter, create a non-background mask, and adjust

% contrast using only non-background pixels.

% --- CODED ANSWER ---

% Load image

imOriginal = imread("Section3.png");

% Convert to grayscale if needed

if ndims(imOriginal) == 3

imOriginal = rgb2gray(imOriginal);

end

% Why: many image processing functions expect grayscale for simple intensity operations

% Display original

figure; imshow(imOriginal);

title('Original Image');

% 1) Median filter (good for salt-and-pepper noise)

imMed = medfilt2(imOriginal);

% Why: median filtering removes isolated bright/dark pixels without blurring edges too much

% Important for salt-and-pepper noise specifically

figure;

imshowpair(imOriginal, imMed, "montage", "scaling", "none");

title('Original (left) vs Median Filtered (right)');

% 2) Averaging filter (for comparison)

avgFilt = (1/9) * ones(3,3);

imAvg = filter2(avgFilt, imOriginal);

imAvg = uint8(imAvg);

% Why: averaging smooths by local mean, but it can blur edges/details

% Important to compare filter behavior and choose the right tool

figure;

imshowpair(imOriginal, imAvg, "montage", "scaling", "none");

title('Original (left) vs Averaging Filtered (right)');

% 3) Histogram of median-filtered image

figure;

imhist(imMed);

title('Histogram of Median-Filtered Image');

% Why: histogram shows how pixel intensities are distributed

% Important for deciding if contrast needs stretching

% 4) Create mask for non-background pixels (background assumed black = 0)

mask = (imMed ~= 0);

% Why: excludes black background so contrast adjustment is based on the tissue/image region

% Important because background can dominate histogram and reduce contrast stretch quality

% 5) Compute contrast limits using only non-background pixels

p = imMed(mask);

inLow = double(min(p))/255;

inHigh = double(max(p))/255;

% Why: get intensity limits from only useful pixels, then normalize to [0,1] for imadjust

% 6) Adjust contrast

imAdj = imadjust(imMed, [inLow inHigh], []);

% Why: stretches intensities to use more display range, improving visibility

% Important for clearer MRI detail

% 7) Histogram after adjustment

figure;

imhist(imAdj);

title('Histogram After Contrast Adjustment');

% 8) Final comparison

figure;

imshowpair(imOriginal, imAdj, "montage", "scaling", "none");

title('Original (left) vs Median + Contrast Adjusted (right)');

%% =========================================================

%% 3.2 Multichannel statistical analysis

%% =========================================================

% PRACTICE QUESTION:

% Given a CSV table with 10 subject columns of heart rate data (plus a first

% column that is time/index):

% 1) Compute mean, median, range, std for each subject

% 2) Find which subject has max and min mean HR

% 3) Compute hours each subject spent above (mean + 20 bpm)

% if each row = 5 seconds

% --- CODED ANSWER ---

% Load and convert table

T = readtable("Section2.csv");

A = table2array(T);

% Why: readtable is good for CSVs; table2array makes math/stats easier column-wise

% Select subject columns (assuming columns 2:11 are subjects)

X = A(:,2:11);

% Why: isolates just the 10 heart-rate signals (one column per subject)

% 1) Column-wise statistics (one value per subject)

mean_data = mean(X, 'omitnan');

median_data = median(X, 'omitnan');

range_data = max(X, [], 1) - min(X, [], 1);

std_data = std(X, 0, 'omitnan');

% Why: summarizes each subject's heart-rate behavior (central tendency + spread)

% Important for comparing subjects objectively

% 2) Find max/min mean HR and subject index

[max_mean, subj_max_mean] = max(mean_data);

[min_mean, subj_min_mean] = min(mean_data);

% Why: max/min gives both the value and which subject had it

% Important because lab asks for subject ID via coding (not by eye)

fprintf('Max mean HR: Subject %d (%.2f bpm)\n', subj_max_mean, max_mean);

fprintf('Min mean HR: Subject %d (%.2f bpm)\n', subj_min_mean, min_mean);

% 3) Time above (mean + 20 bpm), row spacing = 5 s

dt = 5;

hours_over = zeros(1,10);

for s = 1:10

threshold = mean_data(s) + 20;

count_above = sum(X(:,s) > threshold);

% Why: sum(logical) counts how many samples exceed the threshold

hours_over(s) = (count_above * dt) / 3600;

% Why: convert sample count -> seconds -> hours

% Important because lab asks for duration in hours, not sample count

end

% Display all subject results in a table

Subject = (1:10)';

resultsTbl = table(Subject, mean_data', median_data', range_data', std_data', hours_over', ...

'VariableNames', {'Subject','MeanHR','MedianHR','RangeHR','StdHR','HoursAboveMeanPlus20'});

% Why: makes results easy to read/check quickly

disp(resultsTbl);

[Image Analyst]

Yep - that's some code. But so what? Are we supposed to do anything about it? Do you have a specific question about it?

This looks like a homework problem. If you have any questions ask your instructor or read the link below to get started:

How do I get help on homework questions on MATLAB Answers? - MATLAB Answers - MATLAB Central

Obviously, if it is homework, we can't give you the full solution because you're not allowed to turn in our code as your own.

Another resource that MathWorks Central offers is the AI Chat Playground (on the blue banner above). You can paste your problem text into there and get code. However to be ethical you need to make sure that your professor is willing to accept AI generated code as your own. If he/she is against that and wants you to create your own code, and you submit code from AI or other people, you could be running a risk.

Answer 1

[Image Analyst]

Right before

%% Compute PSDs for original signals

[f1, PSD1] = myPSD(sig1, Fs);

[f2, PSD2] = myPSD(sig2, Fs);

you're going to have to come up with a sig1, sig2, and Fs and then run the program. If you want you can just use rand() or a rand vector added to a sine wave vector to give a noisy sinne wave signal. Like

numSamplePoints = 1000;
period = 500;
signalAmplitude = 10;
noiseAmplitude = 3;
maxXValue = 2000;
x = linspace(1, maxXValue, numSamplePoints);
sig1 = signalAmplitude * sin(2 * pi * x / period) + noiseAmplitude * rand(1, numSamplePoints);
plot(x, sig1, 'b-');
grid on;
yline(0, 'LineWidth', 1);
% Now create signal 2 with a different amplitude and phase
signalAmplitude = 14;
sig2 = signalAmplitude * sin(2 * pi * (x - 40) / period) + noiseAmplitude * rand(1, numSamplePoints);
hold on;
plot(x, sig2, 'r-');
legend('sig1', 'y=0', 'sig2')

●

Experiment around with different Fs and other parameters and observe how the fft spectrum changes.