From 14ca4cc8455dbe6eea74cfca2acd9111af3dbb4b Mon Sep 17 00:00:00 2001 From: Guillaume Le chartier Date: Mon, 18 Mar 2024 09:27:00 +0100 Subject: [PATCH] added files --- frequencySpectrum.m | 61 +++++++++++++++++++++++++++++++++++++++++++++ spectrogram.m | 38 ++++++++++++++++++++++++++++ unknownsignal.csv | 1 + 3 files changed, 100 insertions(+) create mode 100644 frequencySpectrum.m create mode 100644 spectrogram.m create mode 100644 unknownsignal.csv diff --git a/frequencySpectrum.m b/frequencySpectrum.m new file mode 100644 index 0000000..8576cb5 --- /dev/null +++ b/frequencySpectrum.m @@ -0,0 +1,61 @@ +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) +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/spectrogram.m b/spectrogram.m new file mode 100644 index 0000000..e26c2c0 --- /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, 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_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*samplingFreq/2/samplingFreq,:)); % magnitude in range 0