Day 3

% Signal parameters
Fs  = 1000;           % sampling rate (Hz)
T   = 1/Fs;           % sample period
N   = 1024;           % number of samples
t   = (0:N-1) * T;    % time vector

% Composite signal: 50 Hz + 120 Hz + noise
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t) + 0.5*randn(1,N);

% Compute FFT
X    = fft(x);
X    = X(1:N/2+1);    % keep positive frequencies only
mag  = abs(X) / N;    % magnitude spectrum
mag(2:end-1) = 2*mag(2:end-1);  % double (single-sided)

% Frequency axis
f = Fs * (0:N/2) / N;

% Plot time domain and frequency domain
figure;
subplot(2,1,1);
plot(t(1:200), x(1:200));    % first 200 ms
title('Time Domain Signal'); xlabel('Time (s)'); ylabel('Amplitude');
grid on;

subplot(2,1,2);
plot(f, mag, 'b', 'LineWidth', 1.5);
title('Single-Sided Amplitude Spectrum');
xlabel('Frequency (Hz)'); ylabel('Amplitude');
xlim([0, 200]);
grid on;

% Mark the peaks
[pks, locs] = findpeaks(mag, f, 'MinPeakHeight', 0.3);
hold on;
scatter(locs, pks, 80, 'r', 'filled');
text(locs+2, pks, string(round(locs)) + " Hz", 'FontSize', 9);

Fs = 1000;   % sampling rate

% ── FIR Low-pass filter (window method) ───────────────────────────
% fir1(order, cutoff_normalized)  — cutoff is fraction of Nyquist (Fs/2)
cutoff_hz = 100;
cutoff_norm = cutoff_hz / (Fs/2);   % normalize to [0,1]
b_fir = fir1(64, cutoff_norm, 'low');  % 64th order FIR

% Frequency response
[H, f] = freqz(b_fir, 1, 1024, Fs);
figure;
subplot(2,1,1);
plot(f, 20*log10(abs(H)));
title('FIR Low-pass Filter'); xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)');
xlim([0, Fs/2]); ylim([-80, 5]); grid on;
xline(cutoff_hz, 'r--', '100 Hz cutoff');

% ── IIR Butterworth filter ─────────────────────────────────────────
% butter(order, cutoff_normalized)
[b_iir, a_iir] = butter(5, cutoff_norm, 'low');
[H2, ~] = freqz(b_iir, a_iir, 1024, Fs);
subplot(2,1,2);
plot(f, 20*log10(abs(H2)));
title('5th Order Butterworth Low-pass'); xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)');
xlim([0, Fs/2]); ylim([-80, 5]); grid on;
xline(cutoff_hz, 'r--', '100 Hz cutoff');

% ── Apply filters to a noisy signal ───────────────────────────────
t_sig = 0:1/Fs:1-1/Fs;
clean = sin(2*pi*50*t_sig);                         % 50 Hz target
noise = 0.8*sin(2*pi*300*t_sig) + 0.3*randn(size(t_sig));  % HF noise
noisy = clean + noise;

filt_fir = filter(b_fir, 1,        noisy);   % FIR filter
filt_iir = filtfilt(b_iir, a_iir, noisy);    % IIR zero-phase

figure;
plot(t_sig(1:100), noisy(1:100), 'k:'); hold on;
plot(t_sig(1:100), filt_fir(1:100), 'b');
plot(t_sig(1:100), filt_iir(1:100), 'r');
legend({'Noisy', 'FIR filtered', 'IIR zero-phase'});
title('Filtering Comparison'); xlabel('Time (s)'); grid on;

Fs = 4000;
t  = 0:1/Fs:2;

% Chirp signal: frequency sweeps from 100 Hz to 1000 Hz
x = chirp(t, 100, 2, 1000);

% Add a burst of 600 Hz from t=0.8 to t=1.2
burst_mask = (t >= 0.8) & (t <= 1.2);
x = x + 0.5 * sin(2*pi*600*t) .* burst_mask;

% Spectrogram using STFT
% spectrogram(x, window, overlap, nfft, Fs)
window  = hann(256);
overlap = 200;
nfft    = 512;

figure;
spectrogram(x, window, overlap, nfft, Fs, 'yaxis');
colormap('jet');
title('Spectrogram: Chirp + Burst');
xlabel('Time (s)'); ylabel('Frequency (Hz)');
colorbar;
ylim([0, 1500]);

% Manual STFT if you need the data (not just the plot)
[S, F, T] = spectrogram(x, window, overlap, nfft, Fs);
power_dB = 10*log10(abs(S).^2 + eps);
fprintf('Spectrogram: %d freq bins × %d time frames\n', size(S,1), size(S,2));