Merge branch 'develop'

frequencySpectrum.m

@@ -0,0 +1,48 @@
+function power = frequencySpectrum(signal, fs)
+%%%%%%%%%%%%%%%%%%
+%function frequencySpectrum(signal, fs)
+%
+% Task: Display the power spectrum of a given signal
+%
+% Input:
+%	- signal: the input signal to process
+%	- fs: the sampling rate
+%
+% Output: 
+%	- power: power spectrum of the signal
+%	
+%
+% Guillaume Gibert, guillaume.gibert@ecam.fr
+% 25/04/2022
+%%%%%%%%%%%%%%%%%%
+
+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)')
+

main.m

@@ -7,17 +7,75 @@
 %Last Modified: 20/04/23 10:08
 %
 %%%%%%%%%%%%%%%%
+pkg load signal;
 
-Fs = 300; %Hz
+samplingFreq = 300; %Hz
+
+fMin = 30; %Hz
+fMax = 40; %Hz
 
 signal = csvread('unknownsignal.csv');
 
-signalDuration = size(signal)/Fs; %s
+signalDuration = size(signal,2)/samplingFreq; %s
-t=[0:1/Fs:(size(signal,2)-1)/Fs];
+t=[0:1/samplingFreq:(size(signal,2)-1)/samplingFreq];
+windowDuration = signalDuration/2;
 
 figure;
 plot(t,signal)
 xlabel('Time (s)');
 ylabel('Sound (dB)');
+title('unknown Signal');
+
+blackmanWin = zeros(1, length(t));
+for l_sample=1:windowDuration*samplingFreq
+	blackmanWin(l_sample+signalDuration*samplingFreq/4) = (0.42 - 0.5 * cos(2*pi*(l_sample)/(signalDuration*samplingFreq/2)) + 0/08*cos(4*pi*(l_sample)/(windowDuration*samplingFreq/2)));
+end
+
+% plot Blackman window
+%~ figure;
+%~ plot(t, blackmanWin);
+
+% apply the Blackman window
+for l_sample=1:length(t)
+	signal_blackman(l_sample) = signal(l_sample) * blackmanWin(l_sample);
+end
+
+% plot signal windowed by rectangular window
+%~ figure;
+%~ plot(t, signal_blackman);
+
+% plot the frequency spectrum of this windowed signal
+power_blackman = frequencySpectrum(signal_blackman, samplingFreq);
 
 
+fft_signal = fft(signal);
+size(fft_signal)
+
+figure;
+plot(t, fft_signal);
+xlabel('Freq (Hz)');
+ylabel('Sound (dB)');
+title('FFT of signal');
+
+
+N=size(signal)(2);
+freq=(0:N-1)*samplingFreq/N;
+%frequency from frame rate
+
+
+minfreq=30; %Hz
+maxfreq=40; %Hz
+idx_min = find(freq >= minfreq, 1);
+idx_max = find(freq <= maxfreq, 1, 'last');
+filtered_freq = freq(idx_min:idx_max);
+figure;
+plot(t, filtered_freq);
+
+
+
+%[a, b] = butter(1, [1/fMax, 1/fMin], 'bandpass');
+%filtered_signal = filter(b,1, signal);
+%figure;
+%plot(t,filtered_signal)
+%xlabel('Time (s)');
+%ylabel('Sound (dB)');