diff --git a/blackmanWin.m b/blackmanWin.m new file mode 100644 index 0000000..1aeb2af --- /dev/null +++ b/blackmanWin.m @@ -0,0 +1,42 @@ +function signal_win = blackmanWin(signal, signal_duration, sampling_freq, pad) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%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 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% create the temporal array +t=-signal_duration/2:1/sampling_freq:(signal_duration/2)-1/sampling_freq; +% window duration is half of signal duration +windowDuration = signal_duration/2; + + +blackmanWin = zeros(1, length(t)); + +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; +title("Blackman Window"); hold on; +plot(t, blackmanWin); + +% apply the Blackman window +for l_sample=1:length(signal) + signal_win(l_sample) = signal(l_sample) * blackmanWin(l_sample); +end + +figure; +title("Original and Windowed signals"); hold on; +plot(signal); hold on; +plot(signal_win); + +frequencySpectrum(signal_win, sampling_freq/2,pad); \ No newline at end of file diff --git a/frequencySpectrum.m b/frequencySpectrum.m new file mode 100644 index 0000000..67548c4 --- /dev/null +++ b/frequencySpectrum.m @@ -0,0 +1,55 @@ +function power = frequencySpectrum(signal, fs, pad) +%%%%%%%%%%%%%%%%%% +%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 +% - pad: pad the signal with zeros to the next power of 2 +% +% Output: +% - power: the power spectrum +% +% +% Guillaume Gibert, guillaume.gibert@ecam.fr +% 25/04/2022 +%%%%%%%%%%%%%%%%%% + +n = length(signal); % number of samples + +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/hammingWin.m b/hammingWin.m new file mode 100644 index 0000000..19bbcae --- /dev/null +++ b/hammingWin.m @@ -0,0 +1,48 @@ +function signal_win = hammingWin(signal, signal_duration, sampling_freq, pad) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%function signal_win = hammingWin(signal, signal_duration, sampling_freq, pad) +% ex.: signal_win = hammingWin(signal, 2, 300, 1) +% +% Inputs: +% - signal: location of the signal to window +% - signal_duration: duration of the signal in seconds +% - sampling_freq: sampling frequency in Hz +% - pad: wether or not add zero padding (0 false; 1 true) +% +% Output: +% - signal_win: signal of interest on which a hamming window was applied +% +% Author: Guillaume Gibert, guillaume.gibert@ecam.fr +% Date: 04/03/2024 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% create the temporal array +t=-signal_duration/2:1/sampling_freq:(signal_duration/2)-1/sampling_freq; +% window duration is half of signal duration +windowDuration = signal_duration/2; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% creates the Hamming time window +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +hammingWin = zeros(1, length(t)); + +for l_sample=1:windowDuration*sampling_freq + hammingWin(l_sample+signal_duration*sampling_freq/4) = (0.54 - 0.46*cos(2*pi*(l_sample)/(signal_duration*sampling_freq/2))); +end + +figure; +title("Hamming Window"); hold on; +plot(t,hammingWin); + +% apply the window on input signal +for l_sample=1:length(t) + signal_win(l_sample) = signal(l_sample) * hammingWin(l_sample); +end + +figure; +title("Original and Windowed signals"); hold on; +plot(t, signal); hold on; +plot(t, signal_win); + +frequencySpectrum(signal_win, sampling_freq/2,pad); \ No newline at end of file diff --git a/hanningWin.m b/hanningWin.m new file mode 100644 index 0000000..d6ed0ef --- /dev/null +++ b/hanningWin.m @@ -0,0 +1,48 @@ +function signal_win = hanningWin(signal, signal_duration, sampling_freq, pad) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%function signal_win = hanningWin(signal, signal_duration, sampling_freq, pad) +% ex.: signal_win = hanningWin(signal, 2, 300, 1) +% +% Inputs: +% - signal: location of the signal to window +% - signal_duration: duration of the signal in seconds +% - sampling_freq: sampling frequency in Hz +% - pad: wether or not add zero padding (0=false; 1=true) +% +% Output: +% - signal_win: signal of interest on which a hanning window was applied +% +% Author: Guillaume Gibert, guillaume.gibert@ecam.fr +% Date: 04/03/2024 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% create the temporal array +t=-signal_duration/2:1/sampling_freq:(signal_duration/2)-1/sampling_freq; +% window duration is half of signal duration +windowDuration = signal_duration/2; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% creates the Hanning time window +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +hanningWin = zeros(1, length(t)); + +for l_sample=1:windowDuration*sampling_freq + hanningWin(l_sample+signal_duration*sampling_freq/4) = (0.5 - 0.5*cos(2*pi*(l_sample)/(signal_duration*sampling_freq/2))); +end + +figure; +title("Hanning Window"); hold on; +plot(t,hanningWin); + +% apply the window on input signal +for l_sample=1:length(t) + signal_win(l_sample) = signal(l_sample) * hanningWin(l_sample); +end + +figure; +title("Original and Windowed signals"); hold on; +plot(t, signal); hold on; +plot(t, signal_win); + +frequencySpectrum(signal_win, sampling_freq/2,pad); \ No newline at end of file diff --git a/rectWin.m b/rectWin.m new file mode 100644 index 0000000..2c83c6d --- /dev/null +++ b/rectWin.m @@ -0,0 +1,48 @@ +function signal_win = rectWin(signal, signal_duration, sampling_freq, pad) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%function signal_win = rectWin(signal, signal_duration, sampling_freq, pad) +% ex.: signal_win = rectWin(signal, 2, 300, 1) +% +% Inputs: +% - signal: location of the signal to window +% - signal_duration: duration of the signal in seconds +% - sampling_freq: sampling frequency in Hz +% - pad: wether or not add zero padding (0=false; 1=true) +% +% Output: +% - signal_win: signal of interest on which a rectangular window was applied +% +% Author: Guillaume Gibert, guillaume.gibert@ecam.fr +% Date: 04/03/2024 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +signal_duration = length(signal)/sampling_freq; +% create the temporal array +t=-signal_duration/2:1/sampling_freq:(signal_duration/2)-1/sampling_freq; +% window duration is half of signal duration +windowDuration = signal_duration/2; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% creates the Rectangular time window +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +rectangularWin = zeros(1, length(t)); + +for l_sample=1:windowDuration*sampling_freq + rectangularWin(l_sample + signal_duration) = 1; +end + +figure; +title("Rectangular Window"); hold on; +plot(t,rectangularWin); + +% apply the rectangular window +for l_sample=1:length(t) + signal_win(l_sample) = signal(l_sample) * rectangularWin(l_sample); +end + +figure; +title("Original and Windowed signals"); hold on; +plot(t, signal); hold on; +plot(t, signal_win); + +frequencySpectrum(signal_win, sampling_freq/2,pad); diff --git a/spectrogram.m b/spectrogram.m new file mode 100644 index 0000000..ca44b95 --- /dev/null +++ b/spectrogram.m @@ -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_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); +xlabel('Time (s)'); +ylabel('Frequency (Hz)');