sebastien.dubois
/
Midterm_Sebastien_Dubois_18-03-2024
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'Develop'

This commit is contained in:
Sébastien DUBOIS 2024-03-18 09:08:55 +01:00
parent 9e921d107e f064e56999
commit 96f268022b
7 changed files with 284 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
README.md
View File

@ -1 +1,6 @@
This is the read me file We have 5 codes here.
The one we need to use is the main one.
To use it we just run it.
It will plot 10 graphs representing the signals, the linear frequency and the log frequency

32
blackmanWin.m Normal file
View File

@ -0,0 +1,32 @@
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(2pi(l_sample)/length(signal)) + 0/08cos(4pi(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);

63
frequencySpectrum.m Normal file
View File

@ -0,0 +1,63 @@
function power = frequencySpectrum(signal, fs, resolution, part_signal)
%%%%%%%%%%%%%%%%%%
%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)
xticks(0:0.1fs:nfs);
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('Plot for signal number', part_signal)
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)')

84
main.m Normal file
View File

@ -0,0 +1,84 @@
%%%%%%%%%%%%%%%%%%%%%%
% 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 = 200; % -> 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);
xticks(0:0.2*fps:length(signal)*fps);
xticklabels(0:0.2:length(signal)/fps);
xlabel('Time (s)');
ylabel('Amplitude (a.u.)');
title('Temporal plot');
% 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 0.5 Hz resolution
frequencySpectrum(signal_1, fps, 0.5, 1);
% 2nd part of the signal with 1 Hz resolution
frequencySpectrum(signal_2, fps, 0.5, 2);
%%%%%SPECTRAL ANALYSIS (RECTANGULAR WINDOW)%%%%%
%plots power spectrum with rectangular window
% 1st part of the signal with 0.5 Hz resolution
zoomFrequencySpectrum(signal_1, fps, 0.5, 1);
% 2nd part of the signal with 1 Hz resolution
zoomFrequencySpectrum(signal_2, fps, 0.5, 2);
%%%%%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, 0.5, 1);
signal_2_win = blackmanWin(signal_2);
% 2nd part of the signal with 1 Hz resolution
frequencySpectrum(signal_2_win, fps, 0.5, 2);
%%%%%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
zoomFrequencySpectrum(signal_1_win, fps, 0.5, 1);
signal_2_win = blackmanWin(signal_2);
% 2nd part of the signal with 1 Hz resolution
zoomFrequencySpectrum(signal_2_win, fps, 0.5, 2);

35
spectrogram.m Normal file
View File

@ -0,0 +1,35 @@
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');
xlim([0 length(signal)/samplingFreq-1/samplingFreq]);
ylabel('Amplitude (norm. unit)');
subplot(2,1,2);
step = fix(step_sizesamplingFreq/1000); % one spectral slice every step_size ms
window = fix(window_sizesamplingFreq/1000); % window_size ms data window
[S, f, t] = specgram(signal);
specgram(signal, 2^nextpow2(window), samplingFreq, window, window-step);
xlabel('Time (s)');
ylabel('Frequency (Hz)');

1
unknownsignal.csv Normal file
View File

File diff suppressed because one or more lines are too long

63
zoomFrequencySpectrum.m Normal file
View File

@ -0,0 +1,63 @@
unction power = zoomFrequencySpectrum(signal, fs, resolution, part_signal)
%%%%%%%%%%%%%%%%%%
%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
%
%
% Sébastien Dubois, sebastien.dubois@ecam.fr
% 18/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)
xticks(0:0.1fs:nfs);
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');
xlim([5, 20]);
title('Zoomed plot for signal number', part_signal)
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)')