Skip to content

Commit 9091232

Browse files
leabouffautleabouffaut
authored
First upload of the new files
0 parents  commit 9091232

6 files changed

+515
-0
lines changed

‎generate_observation_SNR_controlled.m

+85
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
% x = generate_observation_SNR_controlled(noise,signal,fs,fsig,time_of_arrival,SNR)
2+
%
3+
% This fonction can be used to generate an observation with a controlled
4+
% SNR, following this definition SNR = 10*log10(signal_power/noise_power)
5+
% with x_power = x_energy/length(x). This fonction is build for tonals.
6+
% Therefore, the energy term is relative to the energy in the signal's
7+
% frequency band, define by fsig = [fmin fmax]. The noise energy is
8+
% measured at each time_of_arrival value for the duration of the signal,
9+
% and the signal's amplitude is ajusted by the variable coef to satisfy the
10+
% imposed SNR. This function is based on the method proposed (Mellinger &
11+
% Clark 2006, see Mobysound website: http://www.mobysound.org/software.html)
12+
%
13+
% INPUTS:
14+
% - noise, noise vector where the signal is going to be injected
15+
% with a controlled SNR, sampled at fs;
16+
% - signal, the signal that is going to be injected at each time_of_arrival
17+
% sampled at fs;
18+
% - fs, the sampling frequency (Hz);
19+
% - fsig, fsig = [fmin fmax], the signal's frequency band;
20+
% - time_of_arrival, time vector containig the starting time where
21+
% signals will be injected;
22+
% - SNR, imposed Signal to Noise Ratio (dB);
23+
%
24+
% OUTPUT:
25+
% - x, temporal observation vector of a mixture of signal and noise at
26+
% the required SNR
27+
% - real_SNR (optional), double-checked SNR value according
28+
29+
30+
function [x,real_SNR] = generate_observation_SNR_controlled(noise,signal,fs,fsig,time_of_arrival,SNR)
31+
% pre-processing of the data
32+
noise = noise - mean(noise); % so the noise is centered
33+
signal = signal - mean(signal); signal = signal/max(abs(signal)); % the signal is normalized
34+
35+
%% Signal and noise mixture
36+
% Initialize the received signal with the noise
37+
received_signal = noise;
38+
39+
% We want to measure the exact noise power at the TOA where we're going to
40+
% inject the signal, to adjust its power and have the required SNR.
41+
N = length(signal);
42+
N_2 = 2 .^ nextpow2(N);
43+
44+
% Frequency indexes corresponding to fmin and fmax of the simulated Z-call
45+
% +/- 1 Hz (upper and lower)
46+
i0 = round((fsig(1)-1)/(fs/2)*N_2/2) + 1; % index for lower freq bound
47+
i1 = round((fsig(2)+1)/(fs/2)*N_2/2) + 1; % ...and upper
48+
49+
% Signal Energy
50+
Ref_Signal_energy = energy_measurement(signal,[i0 i1],[N N_2]);
51+
52+
% Initialize
53+
real_SNR = [];
54+
% For each emission
55+
for emmission_number = 1: length(time_of_arrival)
56+
57+
% Sample where we'll start injecting the call
58+
arrival_sample = round(time_of_arrival(emmission_number)*fs);
59+
% Noise cut from the arrival sample and for the duration of the call
60+
noise_short = received_signal(arrival_sample:arrival_sample+N-1);
61+
% Calculation of noise_short power in the bandwidth of the Z-call
62+
Noise_energy = energy_measurement(noise_short,[i0 i1],[N N_2]);
63+
% Calculation of the coefficient that will be used to set signal at
64+
% the right amplitude
65+
Signal_wanted_energy = Noise_energy * 10^(SNR/10) - Noise_energy;
66+
coef = sqrt(Signal_wanted_energy/Ref_Signal_energy);
67+
68+
% Received signal
69+
received_signal(arrival_sample:arrival_sample+N-1) = received_signal(arrival_sample:arrival_sample+N-1) + coef*signal;
70+
71+
if nargout == 2
72+
% SNR (to verify)
73+
Observation_energy = energy_measurement(received_signal(arrival_sample:arrival_sample+N-1),[i0 i1],[N N_2]);
74+
real_SNR = [real_SNR 10*log10(Observation_energy/Noise_energy)];
75+
end
76+
77+
end
78+
% Returned built simulated signal
79+
x = received_signal;
80+
x = x-mean(x); x = x/max(x);
81+
82+
function energy = energy_measurement(x,I,NN)
83+
pad = [x zeros(1, NN(2)-NN(1))]; % zero-pad to power-of-2 length
84+
FT = abs(fft(pad));
85+
energy = (1/NN(2)) * sum(FT(I(1):I(2)).^2)/NN(1);

‎lea_Baumgartner.m

+83
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
% function [f0_baumgartner, time_baumgartner, ampl_baumgartner] = lea_Baumgartner(p,f,t,Ts)
2+
%
3+
% This pitch tracking function is based on the algoritm described in
4+
% Baumgartner2011. First a "thresholding is applied. Then, to find the
5+
% pitch tracks, then a cost calculation is realized to associates STFT
6+
% pixels of a track together and have only one f0 value/time.
7+
%
8+
% INPUTS:
9+
% - spectrogram p,f,t
10+
% - Ts: fix threshold in dB.
11+
%
12+
% OUTPUTS:
13+
% - f0_inst_freq: Pitch track from instantaneous frequency estimate,
14+
% - time_inst_freq: time of the HPS pitch track,
15+
% - ampl_inst_freq: amplitude of the tracks.
16+
17+
function [f0_baumgartner, time_baumgartner, ampl_baumgartner] = lea_Baumgartner(p,f,t,Ts)
18+
19+
threshold = 10.^(Ts/10); % real amplitude threshold
20+
21+
% Application of the threshold to the data
22+
p_threshold = NaN(size(p));
23+
for i = 1:length(f)
24+
ind = find(p(i,:)>=threshold);
25+
p_threshold(i,ind) = p(i,ind);
26+
end
27+
28+
% ---------------------- COST CALCULATION -------------------- %
29+
% The cost is calculated in two steps. First for the first "detected" value.
30+
% Then for the following ones. % i is the time index, that starting at the first detected point (i0).
31+
% <=> increase columns.
32+
33+
% Create vectors of indexes for both time and frequency index for
34+
% "detected" calls
35+
[line,col] = find(isnan(p_threshold)==0); % line (= freq = j) and column (= time = i) index of non-NaN values
36+
weight = 20 ; % User define weight
37+
38+
% Empty matrix declaration
39+
Cost = zeros(1,length(t));
40+
Detected_freq = NaN(1,length(t));
41+
Detected_amp = NaN(1,length(t));
42+
% For the first detected value
43+
% i = i0 + 1
44+
i = 2; j = 2; last_j = j;
45+
P = weight * abs(log2(f(line(j-1))/f(line(j))));
46+
Cost(i) = P - p_threshold(line(j),col(i));
47+
48+
% For the rest of the temporal observation
49+
% i > i0 + 1 We go foward in time
50+
for i = 3:length(t)-2
51+
% Are they any "detected frequency" ?
52+
frequency_index = find(isnan(p_threshold(:,i))==0);
53+
% If they are, the following developpment aims to find the frequency
54+
% that minimize the cost, (especially in the case that there are
55+
% several simultaneous detected frequencies).
56+
57+
if isempty(frequency_index) == 0
58+
K = length(frequency_index); % Number of "Detected frequencies"
59+
P = zeros(1,K); % P vector
60+
A = p_threshold(frequency_index,i)'; % Amplitudes
61+
for k = 1:K
62+
P(k) = weight * abs(log2(f(last_j)/f(frequency_index(k))));
63+
end
64+
Cost(i) = min((Cost(i-1) + P - A)); % minimum cost value
65+
if isnan(Cost(i))== 0
66+
ind_freq_Cost_min = find( (Cost(i-1) + P - A) == Cost(i),1); % Corresponding frequency
67+
Detected_freq(i) = f(frequency_index(ind_freq_Cost_min));
68+
Detected_amp(i) = A(ind_freq_Cost_min);
69+
last_j = ind_freq_Cost_min;
70+
end
71+
end
72+
73+
end
74+
75+
for i = length(t)-2:-1:1
76+
end
77+
78+
% % ---------------------- -------------------- -------------------- %
79+
80+
81+
f0_baumgartner = Detected_freq;
82+
time_baumgartner = t;
83+
ampl_baumgartner = Detected_amp;

‎lea_HPS.m

+47
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,47 @@
1+
% function [f0_HPS, time_HPS, ampl_HPS] = lea_HPS(p,f,t)
2+
%
3+
% This pitch tracking function is based on the Harmonic Product Spectrum
4+
% (HPS) algorithm, as described un (Cuadra2001) using the spectrogram as
5+
% time-framing and DFT function.
6+
%
7+
% INPUTS:
8+
% - spectrogram p,f,t
9+
%
10+
% OUTPUTS:
11+
% - f0_HPS: Pitch track from HPS algorithm,
12+
% - time_HPS: time of the HPS pitch track,
13+
% - ampl_HPS: amplitude of the tracks.
14+
15+
function [f0_HPS, time_HPS, ampl_HPS] = lea_HPS(p,f,t)
16+
% function [f0_HPS, time_HPS, ampl_HPS] = lea_HPS(x,fs,fft_size,overlap)
17+
% [~,f,t,p] = spectrogram(x,hann(fft_size),round((overlap/100)*fft_size),fft_size,fs);
18+
19+
p_HPS = ones(length(f),length(t)); % initailization of p_HPS matrix to 1
20+
nb_harmo = 2; % number of times the spctrum is downsampeled. In our case is to low to do it more thant twice.
21+
% For each time bin, the spectrum is downsampeled nb_harmo times. We keep
22+
% the smallest frequency vector size. and complete the rest by nan values
23+
% then the product of p_HPS(:,i) spectrum and the downsampled spectrum is
24+
% realized.
25+
26+
for i = 1:length(t)
27+
h = p(:,i) ;% spectre au temps i
28+
for n = nb_harmo:-1:1
29+
h_down = downsample(h,n);
30+
diff_size = length(f)-length(h_down);
31+
h_down_sized = [h_down; NaN(diff_size,1)];
32+
p_HPS(:,i) = p_HPS(:,i).*h_down_sized;
33+
end
34+
end
35+
[a,~] = size(p_HPS);
36+
f_HPS = linspace(0,50,a)*nb_harmo; %the associated frequency vector
37+
38+
% find Max frequency of the HPS, and the associated amplitude
39+
f0_HPS = zeros(1,length(t));
40+
ampl_HPS = zeros(1,length(t));
41+
42+
for i = 1:length(t)
43+
[M,I] = max(p_HPS(:,i));
44+
f0_HPS(i) = f_HPS(I);
45+
ampl_HPS(i) = sqrt(M);
46+
end
47+
time_HPS = t;

‎lea_inst_freq.m

+56
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,56 @@
1+
% function [f0_inst_freq, time_inst_freq, ampl_inst_freq] = lea_inst_freq(p,f,t,x,fs)
2+
%
3+
% This pitch tracking function is based on instantaneous frequency
4+
% estimate in time. A smoothing median filter is applied as well as
5+
% mean variance measurement for rejecting part of the data.
6+
%
7+
% INPUTS:
8+
% - spectrogram p,f,t
9+
% - x: observation signal,
10+
% - fs: sampling frequency (Hz),
11+
%
12+
% OUTPUTS:
13+
% - f0_inst_freq: Pitch track from instantaneous frequency estimate,
14+
% - time_inst_freq: time of the HPS pitch track,
15+
% - ampl_inst_freq: amplitude of the tracks.
16+
17+
18+
function [f0_inst_freq, time_inst_freq, ampl_inst_freq] = lea_inst_freq(p,f,t,x,fs)
19+
% function [f0_inst_freq, time_inst_freq, ampl_inst_freq] = lea_inst_freq(x,fs,fft_size,overlap)
20+
%[~,f,t,p] = spectrogram(x,hann(fft_size),round((overlap/100)*fft_size),fft_size,fs);
21+
22+
% Instantaneous frequency measurement
23+
z = hilbert(x);
24+
instfreq = fs/(2*pi)*diff(unwrap(angle(z)));
25+
[x] = vector_orientation(x,'line');
26+
27+
% "Smooth" the instantaneous frequency using a median filter
28+
instfreq_filt = medfilt1(instfreq,21);
29+
30+
% Amplitude of these frequencies on the spectrogram
31+
tx = (0:length(x)-1)/fs;
32+
p_interp = interp2(t,f,p,tx,f);
33+
ampl_inst_freq = NaN(size(instfreq_filt));
34+
for i = 1:length(instfreq_filt)
35+
int_f = find(f>= instfreq_filt(i),1);
36+
if isempty(int_f)==0,
37+
ampl_inst_freq(i) = p_interp(int_f,i);
38+
end
39+
end
40+
41+
% win_size = 501; % 5s
42+
% var_f0 = zeros(size(instfreq));
43+
% instfreq_ZP = [zeros(1,floor(win_size/2)) instfreq' zeros(1,floor(win_size/2))];
44+
% for i = 1:length(var_f0),
45+
% temp = instfreq_ZP (i:i+win_size-1);
46+
% var_f0(i) = var(temp(isnan(temp)==0));
47+
% end
48+
%
49+
% f0_inst_freq = instfreq_filt(var_f0<50);
50+
% time_inst_freq = tx(var_f0<50);
51+
% ampl_inst_freq = ampl_inst_freq(var_f0<50);
52+
53+
f0_inst_freq = [instfreq_filt 0];
54+
time_inst_freq = tx;
55+
ampl_inst_freq = [ampl_inst_freq 0];
56+

0 commit comments

Comments
 (0)