77 lines
1.8 KiB
Matlab
77 lines
1.8 KiB
Matlab
## Copyright (C) 2023 Loic
|
|
##
|
|
## This program is free software: you can redistribute it and/or modify
|
|
## it under the terms of the GNU General Public License as published by
|
|
## the Free Software Foundation, either version 3 of the License, or
|
|
## (at your option) any later version.
|
|
##
|
|
## This program is distributed in the hope that it will be useful,
|
|
## but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
## GNU General Public License for more details.
|
|
##
|
|
## You should have received a copy of the GNU General Public License
|
|
## along with this program. If not, see <https://www.gnu.org/licenses/>.
|
|
|
|
## -*- texinfo -*-
|
|
## @deftypefn {} {@var{retval} =} freqRangeSpectrum (@var{input1}, @var{input2})
|
|
##
|
|
## @seealso{}
|
|
## @end deftypefn
|
|
|
|
## Author: Loic <Loic@LAPTOP-1F7SF0F7>
|
|
## Created: 2023-04-13
|
|
|
|
function [val, ind, f] = freqRangeSpectrum (signal, fs, f1, f2, formant, pad)
|
|
n = length(signal);
|
|
if (pad)
|
|
n = 2^nextpow2(n);
|
|
end
|
|
|
|
f = (0:n-1)*(fs/n);
|
|
|
|
idx = find(f >= f1 & f <= f2); %define the index of the freq range
|
|
f = f(idx);
|
|
|
|
tic
|
|
y = fft(signal, n);% compute DFT of input signal
|
|
duration = toc
|
|
|
|
power = abs(y).^2/n;
|
|
power = power(idx);%limit the signal to the frequency ROI
|
|
[val, ind] = max(power);
|
|
|
|
%lowpass for the formant, moving average
|
|
|
|
for j = 1:length(idx)
|
|
tot = 0;
|
|
for k = j-18:j
|
|
if k <=0
|
|
tot = tot + power(1);
|
|
else
|
|
tot = tot + power(k);
|
|
endif
|
|
endfor
|
|
power_avg(j) = 1/18*(tot);
|
|
endfor
|
|
|
|
figure;
|
|
|
|
subplot(1,2,1) % time plot
|
|
if (pad)
|
|
signal = [ signal, zeros( n-length(signal), 1)'];
|
|
end
|
|
|
|
plot(0:1/fs:(n-1)/fs,signal);
|
|
xlabel('Time (s)');
|
|
ylabel('Amplitude (a.u.)');
|
|
|
|
subplot(1,2,2) % freq range plot
|
|
plot(f,10*log10(power/power(ind)));
|
|
if formant == 1
|
|
plot(f, 10*log10(power_avg/power(ind)), 'r');
|
|
endif
|
|
xlabel('Frequency (Hz)')
|
|
ylabel('Power (dB)')
|
|
endfunction
|