From 893829e257d6cd9d7ad705d6e46b962f4e7a0080 Mon Sep 17 00:00:00 2001
From: "ly.pechvattana" <ly.pechvattana@ecam.fr>
Date: Thu, 20 Apr 2023 21:17:14 +0700
Subject: [PATCH] submission

---
 frequencySpectrum.m | 32 ++++++++++++++++++++++++++++++++
 main.m              | 29 +++++++++++++++++++++++++++++
 spectrogram.m       | 38 ++++++++++++++++++++++++++++++++++++++
 3 files changed, 99 insertions(+)
 create mode 100644 frequencySpectrum.m
 create mode 100644 spectrogram.m

diff --git a/frequencySpectrum.m b/frequencySpectrum.m
new file mode 100644
index 0000000..61527f0
--- /dev/null
+++ b/frequencySpectrum.m
@@ -0,0 +1,32 @@
+function power = frequencySpectrum(signal, fs)
+
+n = length(signal);          % number of samples
+
+y = fft(signal, n);% compute DFT of input signal
+power = abs(y).^2/n;    % power of the DFT
+
+[val, ind] = max(power); % find the mx value of DFT and its index
+
+% plots
+figure;
+
+subplot(1,3,1) % time plot 
+t=0:1/fs:(n-1)/fs; % time range
+plot(t, signal)
+xticks(0:0.1*fs:n*fs);
+xticklabels(0:0.1:n/fs);
+xlabel('Time (s)');
+ylabel('Amplitude (a.u.)');
+
+subplot(1,3,2) % linear frequency plot
+f = (0:n-1)*(fs/n);     % frequency range
+plot(f,power, 'b*'); hold on;
+plot(f,power, 'r');
+xlabel('Frequency (Hz)')
+ylabel('Power (a.u.)')
+
+subplot(1,3,3) % log frequency plot
+plot(f,10*log10(power/power(ind)));
+xlabel('Frequency (Hz)')
+ylabel('Power (dB)')
+
diff --git a/main.m b/main.m
index 468f957..82a4865 100644
--- a/main.m
+++ b/main.m
@@ -10,9 +10,38 @@ pkg load signal
 %load data
 signal = csvread('unknownsignal.csv');
 
+%input data
+fs = 1000; %Sampling freq (Hz)
+fc = 100; %Cut-off Freq (Hz)
+n = 4; % Filter order
+
+samplingDuration = length(signal);
+window_size = samplingDuration/2;
+
 %plotting signal in time domain
 figure;
 plot(signal);
 xlabel('Time in [s]');
 ylabel('Amplitude');
 title('Signal data in time domain');
+
+%frequency Spectrum
+frequencySpectrum(signal,samplingFreq);
+
+%spectrogram
+%spectrogram(signal, samplingFreq, 30, 5)
+
+%low pass filter
+[b, a] = butter(n, fc/(fs/2), 'low');
+signal_filtered = filter(b, a, signal);
+t = 0 : length(signal) - 1;
+subplot(2,1,1);
+plot(t, signal);
+title('Original Signal');
+xlabel('Time (samples)');
+ylabel('Amplitude');
+subplot(2,1,2);
+plot(t, signal_filtered);
+title('Low-pass Filtered Signal');
+xlabel('Time (samples)');
+ylabel('Amplitude');
\ No newline at end of file
diff --git a/spectrogram.m b/spectrogram.m
new file mode 100644
index 0000000..2e578ac
--- /dev/null
+++ b/spectrogram.m
@@ -0,0 +1,38 @@
+function spectrogram(signal, samplingFreq, step_size, window_size)
+%%%%%%%%%%%%%%%%%%%%%%%
+%function spectrogram(signal, samplingFreq, step_size, window_size)
+% ex.:  spectrogram(signal, samplingFreq, step_size, window_size)
+%
+% Task: Plot the spectrogram of a given signal
+%
+% Inputs:
+%	-signal: temporal signal to analyse 
+%	-samplingFreq: sampling frequency of the temporal signal
+% 	-step_size: how often the power spectrum will be computed in ms
+%	-window_size: size of the analysing window in ms
+%
+% Ouput: None
+%
+% author: Guillaume Gibert (guillaume.gibert@ecam.fr)
+% date: 14/03/2023
+%%%%%%%%%%%%%%%%%%%%%%%
+
+figure;
+	subplot(2,1,1);
+t=0:1/samplingFreq:length(signal)/samplingFreq-1/samplingFreq;
+plot(t, signal');
+xlim([0 length(signal)/samplingFreq-1/samplingFreq]);
+ylabel('amplitude (norm. unit)');
+	subplot(2,1,2);
+step = fix(step_size*samplingFreq/1000);     % one spectral slice every step_size ms
+window = fix(window_size*samplingFreq/1000);  % window_size ms data window
+fftn = 2^nextpow2(window); % next highest power of 2
+[S, f, t] = specgram(signal, fftn, samplingFreq, window, window-step);
+S = abs(S(2:fftn*4000/samplingFreq,:)); % magnitude in range 0<f<=4000 Hz.
+S = S/max(S(:));           % normalize magnitude so that max is 0 dB.
+S = max(S, 10^(-40/10));   % clip below -40 dB.
+S = min(S, 10^(-3/10));    % clip above -3 dB.
+imagesc (t, f, log(S));    % display in log scale
+set (gca, "ydir", "normal"); % put the 'y' direction in the correct direction
+xlabel('time (s)');
+ylabel('frequency (Hz)');
\ No newline at end of file