From 35cdac9c0a23792981a8032c4b94dbc385fd8eb9 Mon Sep 17 00:00:00 2001 From: Gabri6 Date: Thu, 20 Apr 2023 10:39:51 +0200 Subject: [PATCH 1/4] added blackman windowing --- main.m | 30 +++++++++++++++++++++++++++--- 1 file changed, 27 insertions(+), 3 deletions(-) diff --git a/main.m b/main.m index 8f8fe9e..cc8bc0c 100644 --- a/main.m +++ b/main.m @@ -8,16 +8,40 @@ % %%%%%%%%%%%%%%%% -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)'); +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); From 45a36e3e6b4a099dfeeacff8fa2b5d1e4cb29cf9 Mon Sep 17 00:00:00 2001 From: Gabri6 Date: Thu, 20 Apr 2023 11:00:57 +0200 Subject: [PATCH 2/4] fft of the signal --- main.m | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/main.m b/main.m index cc8bc0c..2860a64 100644 --- a/main.m +++ b/main.m @@ -7,6 +7,7 @@ %Last Modified: 20/04/23 10:08 % %%%%%%%%%%%%%%%% +pkg load signal; samplingFreq = 300; %Hz @@ -23,6 +24,7 @@ figure; plot(t,signal) xlabel('Time (s)'); ylabel('Sound (dB)'); +title('unknown Signal'); blackmanWin = zeros(1, length(t)); for l_sample=1:windowDuration*samplingFreq @@ -45,3 +47,20 @@ end % 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('Time (s)'); +ylabel('Sound (dB)'); +title('FFT of signal'); + + +%b = fir1(600, [1/fMax, 1/fMin], 'bandpass'); +%filtered_signal = filter(b, signal); +%figure; +%plot(t,filtered_signal) +%xlabel('Time (s)'); +%ylabel('Sound (dB)'); From 507f249a7c5462dfac6c4ac0ecdf881485df5a60 Mon Sep 17 00:00:00 2001 From: Gabri6 Date: Thu, 20 Apr 2023 11:14:54 +0200 Subject: [PATCH 3/4] add frequency filter --- main.m | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/main.m b/main.m index 2860a64..c7acb7a 100644 --- a/main.m +++ b/main.m @@ -53,13 +53,28 @@ size(fft_signal) figure; plot(t, fft_signal); -xlabel('Time (s)'); +xlabel('Freq (Hz)'); ylabel('Sound (dB)'); title('FFT of signal'); -%b = fir1(600, [1/fMax, 1/fMin], 'bandpass'); -%filtered_signal = filter(b, 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)'); From fe60a1a27d307f4d00479fe92937c8287a9df880 Mon Sep 17 00:00:00 2001 From: Gabri6 Date: Thu, 20 Apr 2023 11:15:41 +0200 Subject: [PATCH 4/4] add external code to see better the frequencies --- frequencySpectrum.m | 48 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 frequencySpectrum.m 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)') +