diff --git a/frequencySpectrum.m b/frequencySpectrum.m new file mode 100644 index 0000000..94617d5 --- /dev/null +++ b/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)') + diff --git a/main.m b/main.m index 8f8fe9e..c7acb7a 100644 --- a/main.m +++ b/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 -t=[0:1/Fs:(size(signal,2)-1)/Fs]; +signalDuration = size(signal,2)/samplingFreq; %s +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)');