first commit

blackmanWin.m

@@ -0,0 +1,36 @@
+function signal_win = blackmanWin(signal)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%function signal_win = blackmanWin(signal)
+%
+% Inputs:
+%	- signal: signal of interest
+%
+% Output:
+%	- signal_win: signal of interest on which a blackman window was applied
+%
+% Author: Guillaume Gibert, guillaume.gibert@ecam.fr
+% Date: 15/03/2024
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+blackmanWin = zeros(1, length(signal));
+for l_sample=1:length(signal)
+	blackmanWin(l_sample) = (0.42 - 0.5 * cos(2*pi*(l_sample)/length(signal)) + 0/08*cos(4*pi*(l_sample)/length(signal)));
+end
+
+% plot Blackman window
+%~ figure;
+%~ plot(blackmanWin);
+
+% apply the Blackman window
+for l_sample=1:length(signal)
+	signal_win(l_sample) = signal(l_sample) * blackmanWin(l_sample);
+end
+
+%~ figure;
+%~ plot(signal);
+%~ hold on;
+%~ plot(signal_win);
+
+
+
+

frequencySpectrum.m

@@ -0,0 +1,64 @@
+function power = frequencySpectrum(signal, fs, resolution)
+%%%%%%%%%%%%%%%%%%
+%function power = frequencySpectrum(signal, fs, pad)
+%
+% Task: Display the power spectrum (lin and log scale) of a given signal
+%
+% Input:
+%	- signal: the input signal to process
+%	- fs: the sampling rate in Hz
+%	- resolution: frequency resolution in Hz, signal will be padded with zeros if necessary
+%
+% Output: 
+%	- power: the power spectrum
+%	
+%
+% Guillaume Gibert, guillaume.gibert@ecam.fr
+% 15/03/2024
+%%%%%%%%%%%%%%%%%%
+
+n = length(signal);          % number of samples
+current_resolution = fs / n;
+if (resolution < current_resolution)
+	n_original = n;
+	n = fs / resolution;
+	signal = [signal zeros(1, n-n_original)];
+end
+
+%~ if (pad)
+	%~ n_original = n;
+	%~ n = 2^(nextpow2(n));
+	%~ signal = [signal zeros(1, n-n_original)];
+%~ end
+
+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)
+title("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');
+title("Frequency Linear ");
+xlabel('Frequency (Hz)')
+ylabel('Power (a.u.)')
+
+subplot(1,3,3) % log frequency plot
+plot(f,10*log10(power/power(ind)));
+title("Frequency Log");
+xlabel('Frequency (Hz)')
+ylabel('Power (dB)')
+

main.m

@@ -0,0 +1,79 @@
+%%%%%%%%%%%%%%%%%%%%%%
+% UNKNOWN SIGNAL
+% Sampling frequency: 300 Hz
+% Duration; 2 s
+% First second: 0.1Hz, 30 Hz, 30.5 Hz, 60 Hz, 61 Hz
+% Second second: 0.1Hz, 32 Hz, 36 Hz, 64 Hz, 72 Hz
+%%%%%%%%%%%%%%%%%%%%%%
+
+% loads the signal package on Octave
+pkg load signal
+
+% loads signal and its characteristics
+signal = csvread('unknownsignal.csv');
+
+%%%%%SIGNAL CHARACTERISTICS%%%%%
+% sets sampling frequency
+fps = 300; % -> freqMax of the signal should be < 150 Hz (Shannon-Nyquisit theorem), in practice freqMax < 60 Hz would be better
+
+% computes the duration of the signal
+duration = length(signal) / fps; % in s
+
+% estimates its original frequency resolution
+resolution = fps /  length(signal); % in Hz
+
+%%%%%STATIONARITY%%%%%
+% temporal plot
+figure;
+plot(signal);
+title("Raw Signal");
+xticks(0:0.2*fps:length(signal)*fps);
+xticklabels(0:0.2:length(signal)/fps);
+xlabel('Time (s)');
+ylabel('Amplitude (a.u.)');
+
+% spectrogram
+step_size = 50; %ms
+window_size = 100; %ms
+spectrogram(signal, fps, step_size, window_size);
+
+% ccl: signal is not stationary, it is composed of 2 parts
+
+%%%%%SPLIT SIGNAL INTO 2 PARTS%%%%%
+% First part: [0 1s]
+signal_1 = signal(1:end/2);
+% Second part: [1s 2s]
+signal_2 = signal(end/2+1:end);
+
+%%%%%SPECTRAL ANALYSIS (RECTANGULAR WINDOW)%%%%%
+%plots power spectrum with rectangular window
+% 1st part of the signal with 1 Hz resolution
+frequencySpectrum(signal_1, fps, 1);
+% 1st part of the signal with 0.5 Hz resolution
+frequencySpectrum(signal_1, fps, 0.5);
+
+% 2nd part of the signal with 1 Hz resolution
+frequencySpectrum(signal_2, fps, 1);
+% 2nd part of the signal with 0.5 Hz resolution
+frequencySpectrum(signal_2, fps, 0.5);
+
+
+
+%%%%%SPECTRAL ANALYSIS (BLACKMAN WINDOW)%%%%%
+%plots power spectrum with blackman window
+signal_1_win = blackmanWin(signal_1);
+% 1st part of the signal with 1 Hz resolution
+frequencySpectrum(signal_1_win, fps, 1);
+% 1st part of the signal with 0.5 Hz resolution
+frequencySpectrum(signal_1_win, fps, 0.5);
+
+signal_2_win = blackmanWin(signal_2);
+% 2nd part of the signal with 1 Hz resolution
+frequencySpectrum(signal_2_win, fps, 1);
+% 2nd part of the signal with 0.5 Hz resolution
+frequencySpectrum(signal_2_win, fps, 0.5);
+
+
+
+
+

spectrogram.m

@@ -0,0 +1,37 @@
+function spectrogram(signal, samplingFreq, step_size, window_size)
+%%%%%%%%%%%%%%%%%%%%%%%
+%function spectrogram(signal, samplingFreq, step_size, window_size)
+% ex.:  spectrogram(signal, 300, 50, 1000)
+%
+% 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');
+title("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
+
+
+[S, f, t] = specgram(signal);
+specgram(signal, 2^nextpow2(window), samplingFreq, window, window-step);
+title("Spectrogram")
+xlabel('Time (s)');
+ylabel('Frequency (Hz)');