如何找出周期性声音信号的频率？

Question

我正在研究步行模式的声音信号，它有明显的规则模式：

​

​

然后我想我可以用FFT函数得到步行的频率(大约从图像中得到1.7Hz )：

    x = walk_5; % Walking sound with a size of 711680x2 double
    Fs = 48000; % sound frquency
    L=length(x); 

    t=(1:L)/Fs; %time base
    plot(t,x);
    figure;

    NFFT=2^nextpow2(L);      
    X=fft(x,NFFT);       
    Px=X.*conj(X)/(NFFT*L); %Power of each freq components       
    fVals=Fs*(0:NFFT/2-1)/NFFT;      
    plot(fVals,Px(1:NFFT/2),'b','LineSmoothing','on','LineWidth',1);         
    title('One Sided Power Spectral Density');       
    xlabel('Frequency (Hz)')         
    ylabel('PSD');

但这并没有给我我所期望的：

FFT结果：

​

​

变焦图像有很多噪声：

​

在1.7Hz附近没有任何信息

​

下面是日志域中的图，使用

    semilogy(fVals,Px(1:NFFT));

不过，它相当对称：

​

我没发现我的密码有什么问题。您有什么解决方案可以轻松地从步行模式中提取出1.7Hz吗？

下面是mat sound.mat?dl=0中音频文件的链接

非常感谢!

Kai

Answer 1

下面是一个很好的解决方案：

在进行FFT之前，获取原始数据的绝对值。这些数据含有大量的高频噪声，无论信号中存在什么低频周期，都会淹没这些噪声。高频噪声的幅度每1.7秒就会变大，幅度的增加对眼睛来说是可见的，而且是周期性的，但是当你用低频正弦波乘以信号，加上所有的东西，你仍然会得到接近于零的东西。取绝对值改变这一点，使这些振幅调制周期在低频。

尝试下面的代码来比较规则数据的FFT和abs(数据)的FFT。注意，我对您的代码采取了一些随意的做法，例如将我假设的两个立体声通道合并成一个单一频道。

x = (walk_5(:,1)+walk_5(:,2))/2; % Convert from sterio to mono
Fs = 48000; % sampling frquency
L=length(x); % length of sample
fVals=(0:L-1)*(Fs/L); % frequency range for FFT

walk5abs=abs(x); % Take the absolute value of the raw data

Xold=abs(fft(x)); % FFT of the data (abs in Matlab takes complex magnitude)
Xnew=abs(fft(walk5abs-mean(walk5abs))); % FFT of the absolute value of the data, with average value subtracted

figure;
plot(fVals,Xold/max(Xold),'r',fVals,Xnew/max(Xnew),'b')
axis([0 10 0 1])
legend('old method','new method')

[~,maxInd]=max(Xnew); % Index of maximum value of FFT
walkingFrequency=fVals(maxInd) % print max value

并为旧方法和新方法绘制从0到10 Hz的FFT，得到如下结果：

​

​

正如你所看到的，它检测到一个大约1.686赫兹的峰值，对于这些数据来说，这是快速傅立叶变换频谱中的最高峰值。

Answer 2

我建议您忘记DFT方法，因为您的信号不适合这种类型的分析，因为许多原因。即使通过观察你感兴趣的频率范围内的频谱，也没有一种简单的估计峰值的方法：

​

​

当然，你可以尝试与PSD/STFT和其他时髦的方法，但这是一个过头。对于这个任务，我可以想到两种相当简单的方法。

第一种是简单地基于自相关函数。

1. 计算ACF
2. 定义它们之间的最小距离。既然你知道预期频率在1.7Hz左右，那么它就相当于0.58s。让我们把它定为0.5秒作为最小距离。
3. 计算发现的峰值之间的平均距离。

这给了我大约1.72赫兹的频率。

​

​

第二种方法是基于对你的信号的观察，你的信号已经有一些周期性的峰值。因此，我们可以使用findpeaks函数简单地搜索它们。

1. 以与以前相同的方式定义最小峰值距离。
2. 定义最小峰值高度。例如，最大峰值的10%。
3. 得到平均差额。

这给了我1.7赫兹的平均频率。

​

​

简单快捷的方法。显然，有一些东西可以改进，例如：

• 精炼阈值
• 同时寻找正负峰
• 注意一些缺失的峰值，例如，由于低振幅

无论如何，这应该让你开始，而不是被困在糟糕的FFT和懒惰的半规范。

代码片段：

load walk_sound

fs = 48000;
dt = 1/fs;

x = walk_5(:,1);
x = x - mean(x);
N = length(x);
t = 0:dt:(N-1)*dt;

% FFT based
win = hamming(N);
X = abs(fft(x.*win));
X = 2*X(1:N/2+1)/sum(win);
X = 20*log10(X/max(abs(X)));
f = 0:fs/N:fs/2;

subplot(2,1,1)
plot(t, x)
grid on
xlabel('t [s]')
ylabel('A')
title('Time domain signal')

subplot(2,1,2)
plot(f, X)
grid on
xlabel('f [Hz]')
ylabel('A [dB]')
title('Signal Spectrum')

% Autocorrelation
[ac, lag] = xcorr(x);
min_dist = ceil(0.5*fs);
[pks, loc] = findpeaks(ac, 'MinPeakDistance', min_dist);

% Average distance/frequency
avg_dt = mean(gradient(loc))*dt;
avg_f = 1/avg_dt;

figure
plot(lag*dt, ac);
hold on
grid on
plot(lag(loc)*dt, pks, 'xr')
title(sprintf('ACF - Average frequency: %.2f Hz', avg_f))


% Simple peak finding in time domain
[pkst, loct] = findpeaks(x, 'MinPeakDistance', min_dist, ...
                            'MinPeakHeight', 0.1*max(x));

avg_dt2 = mean(gradient(loct))*dt;
avg_f2 = 1/avg_dt2;

figure
plot(t, x)
grid on
hold on
plot(loct*dt, pkst, 'xr')
xlabel('t [s]')
ylabel('A')
title(sprintf('Peak search in time domain - Average frequency: %.2f Hz', avg_f2))