如何通过 MATLAB 生成由不同频率滑行段组成的声音？

Question

我想通过 MATLAB 生成声音，该声音由预定义数量的不同频率滑移（扫描）组成。为此我写了一段MATLAB代码。但是，我遇到了两个问题：

1) 当我播放声音时，声音会在整个播​​放期间发出咔哒声。

-> 这可能与前一个扫描段的结尾和下一个扫描段的开头不同的相位角有关。我试图解决这个问题（请参阅下面的代码） - 到目前为止没有成功。您将在此处看到这种声音的频谱图：

spectrogram of concatenated sound which is perceived with several 'soft' clicks

2) 当我生成声音时，有时除了这些更柔和的咔嗒声之外，还会有更明显的咔嗒声。这在频谱图中清晰可见。 --> 在这里我不确定问题可能是什么以及如何避免它。

spectrogram of concatenated sound with additional distinct click

我如何生成声音的代码如下：

clear all;
close all;

%% define stimulus parameters
soundDuration = 1200; % duration of sound
sf = 44100; % sampling rate
ampl = 0.05; % 0.05; % ampl
segmentDuration = 25; % duration of one standard segment in ms
nSegments = soundDuration/segmentDuration; % number of segments of which the sound should consist of
t = 0:1/sf:(0.025-1/sf); % time vector for segment

%% generate sound consisting of n sweep-segments
complexSound = [];
for iSeg = 1:nSegments

    f1 = 1000:10:3000; 
    f1 = randsample(f1,1); % start freq in Hz for current sweep segment
    f2 = 1500:10:4500;
    f2 = randsample(f2,1); % end freq in Hz  for current sweep segment

    if iSeg == 1
        sweep = ampl * chirp(t,f1,segmentDuration/1000,f2,'logarithmic'); % generate sweep-segment withou considering the phase

    else
        sweep = ampl * chirp(t+1/sf,f1,segmentDuration/1000,f2,'logarithmic',ph); % the current sweep starts with a t+1/sf later and with the phase angle with which the previous sweep ended

    end
    ph = -90+360*(f2*t(end)+1/sf); % calculate the phase at the time point at which the current sweep ends and from that calculate the starting phase for the next sweep

    sweep = sweep';
    complexSound = [complexSound; sweep]; % concatenate sweep segments to form the complex sound

end
stim = complexSound; 
sound(stim,sf);

感谢您对解决这些问题的任何帮助。

Answer 1

所以，基本上你想要的是频率调制。好吧，不是真的，但我们基本上可以滥用matlab为其提供的功能。诀窍是这样的：

1. 我们创建一个向量，其元素与最终信号一样多。向量的每个元素都保存最终信号在该特定样本处应具有的频率。
2. 我们使用 fmmod 函数创建一个 1Hz 信号，该信号由我们在 1 中创建的矢量调制。

我尝试调整您的代码并使其尽可能相似。但是，现在所有持续时间都以秒为单位。

请注意，我假设各个段之间的边缘是可以的。因此，例如，如果您有一个段以 1500Hz 结束，而下一个段以 2000Hz 开始，则这两者之间没有平滑的进展。也许 fmmod 函数可以稍微平滑一下，但我没有检查。

但是，这里是代码：

%% clear
clear all;
close all;

%% define stimulus parameters
soundDuration = 1.2; % duration of sound in seconds
sf = 44100; % sampling rate
ampl = 0.05; % 0.05; % ampl
segmentDuration = 0.025; % duration of one standard segment in s
nSegments = round(soundDuration/segmentDuration); % number of segments of which the sound should consist of
samples_per_segment = floor(sf * segmentDuration);

%% generate sound consisting of n sweep-segments
modulator = zeros(1, nSegments * samples_per_segment);
first_idx = 1;
for iSeg = 1:nSegments

    f1 = 1000:10:3000; 
    f1 = randsample(f1,1); % start freq in Hz for current sweep segment
    f2 = 1500:10:4500;
    f2 = randsample(f2,1); % end freq in Hz  for current sweep segment

    modulator(first_idx:first_idx + samples_per_segment-1) = logspace(log10(f1), log10(f2), samples_per_segment); % we add the logarithmic progression from f1 to f2 here
    first_idx = first_idx + samples_per_segment;

end
%% create final sound
stim = fmmod(modulator, 1, sf, 1); % here we abuse the frequency modulation function of matlab.
                                   % we basically tell it to create a 1Hz
                                   % sinewave and modulate it by the vector
                                   % we created in the previous step. the
                                   % last argument tells it basically that
                                   % if i give it a 1, it should increase
                                   % the frequency by 1Hz, if i give it a
                                   % 100, it should increase the freq by
                                   % 100Hz and so on...
stim = ampl .* stim; % multiply by the amplitude

%% plot...
spectrogram(stim,1000,[],[],sf,'yaxis'); % this provides a nice tf plot....

%% play
sound(stim,sf);

Answer 2

我从 Dik J. Hermes 那里收到了解决此问题的另一个解决方案，我也会分享。这里的优点是不依赖 fmmod 函数（上述解决方案所需），该函数仅在商业工具箱中可用。

原始代码（位于本页顶部）中的问题是阶段的管理不正确。下面的代码首先计算瞬时频率，insFreq(t)。如果 stim(t) = a*sin(phi(t))，则 instFreq = phi'(t)/2*pi，或 phi(t) = 2*pi*integral(instFreq)。 instFreq(t) 在 for 循环中计算。在 for 循环之后，通过积分 instFreq 并乘以 2*pi 来计算 phi(t)。通过这种方式，人们注意信号刺激本身是连续的并且不包含相位跳跃。由于瞬时频率的突然变化，刺激的时间导数在从一个片段到另一个片段的过渡处是不连续的。这会导致一些小的基本点击，这是无法避免的（至少对于这里使用的参数（f1、f2、段持续时间）来说不是。

clear all;
close all;

% define stimulus parameters
soundDuration = 1.2; % duration of sound
sf = 44100; % sampling rate
dt = 1/sf;
ampl = 0.05; % 0.05; % ampl
segmentDuration = 0.025; % duration of one segment in s
nSegments = round(soundDuration/segmentDuration); % number of segments of which the sound should consist of
tSegment = 0:dt:(segmentDuration-dt); % time vector for segment

% generate sound consisting of n sweep-segments
complexSound = [];
logInstFreq = [];
startFreqArray = 1000:10:3000;
endFreqArray = 1500:10:4500;
p = tSegment/tSegment(end);
for iSeg = 1:nSegments
    f1 = randsample(startFreqArray,1); % start freq in Hz for current sweep segment
    f2 = randsample(endFreqArray,1); % end freq in Hz  for current sweep segment
    logf1 = log10(f1);
    logf2 = log10(f2);
    logInstFreq = [logInstFreq (logf2-logf1)*p+logf1]; 
end
instFreq = 10.^logInstFreq;
phi = 2*pi*cumsum(instFreq)*dt;
stim = ampl*sin(phi); 
t = 0:dt:(length(instFreq)-1)*dt;
plot(t, instFreq)
sound(stim,sf);
audiowrite('sound.wav', stim, sf);