如何优化需要多个带通滤波器的基于 matlab 的音频可视化工具?
How can I optimize a matlab based audio visualizer that requires multiple bandpass filters?
问题
我接到一项任务,要在 matlab 中使用 FFT 设计一个音频可视化工具。我已经完成了大部分我想用音频可视化工具做的事情,但我 运行 正在顺利地解决 运行 程序的问题,因为它被我的处理器限制了。我想用我的可视化工具做更多的事情,但我不想再实现任何功能,然后才能很好地掌握 运行 的基础知识。
背景
该程序实际上分为两部分。第一部分是音频播放器和 运行 一个名为 CLICK_ME.m
的 matlab 脚本。第二个是currentfft()
。它从歌曲中获取信息,在我厌倦它之前,它绘制了(频率,功率)。然而,这是最基本的可视化工具,我的教授希望我们创建一个独特的音频可视化工具。
我的想法是绘制其中一个统计五边形,其中五个点中每个点的半径是在五个频段之间过滤的功率。
自然地,我实现了 matlab 的 bandpass()
函数,这实现了它的设计目的。然而,它需要的处理能力比我预期的要多得多,随后它遇到了瓶颈,看起来真的很不稳定。
我尝试通过使用 for 循环每隔一个或每隔三个值来减少必须执行计算的点数。我已经这样做了,直到商数减少了 10 为止,但没有取得明显的成功。如果我删除超过十分之一的点,我会从 matlab 中得到错误,它不能对少于六个样本进行 bandpass()
。
我的教授建议使用 gaussmf()
或 rectangularPulse()
来模拟接近带通滤波器的东西,但我不确定如何实现它们或者它们是否会更快。
function currentfft ( player, Y, FS )
sampleNumber = get( player, 'CurrentSample' );
timerVal = get( player, 'TimerPeriod' );
%Get channel one values for our window around the current sample number
s1 = Y(floor(sampleNumber-((timerVal*FS)/2)):floor(sampleNumber+((timerVal*FS)/2)),1);
n = length(s1);
p = fft(s1); % take the fourier transform
nUniquePts = ceil((n+1)/2);
p = p(1:nUniquePts); % select just the first half since the second half
% is a mirror image of the first
p = abs(p); % take the absolute value, or the magnitude
p = p/n; % scale by the number of points so that
% the magnitude does not depend on the length
% of the signal or on its sampling frequency
p = p.^2; % square it to get the power
% multiply by two
if rem(n, 2) % odd nfft excludes Nyquist point
p(2:end) = p(2:end)*2;
else
p(2:end -1) = p(2:end -1)*2;
end
% reduce the number of points actually being filtered
q = 1;
s = 1;
d = int16( length(p) );
pNew = zeros( [ d/10, 1 ] );
while q < d
pNew(s) = p(q);
q = q + 10;
s = s + 1;
end
% try gaussian or rect functions instead of bandpass?
% radius of each section of the pentagon
p0 = abs( bandpass(pNew, [ 1 60 ], FS) );
p1 = abs( bandpass(pNew, [ 60 250 ], FS) );
p2 = abs( bandpass(pNew, [ 250 2e3 ], FS) );
p3 = abs( bandpass(pNew, [ 2e3 8e3 ], FS) );
p4 = abs( bandpass(pNew, [ 8e3 20e3 ], FS) );
% length( p0 )
% length( p1 )
% length( p2 )
% length( p3 )
% length( p4 )
pArr = [ p0, p1, p2, p3, p4 ];
%freqArray = (0:nUniquePts-1) * (FS / n); % create the frequency array
thetaArr = [ pi/2, 4.5*pi/5, 6.5*pi/5, 8.5*pi/5, 10.5*pi/5 ];
% calculating x and y from the radii
x = cos(thetaArr) ./ pArr;
y = sin(thetaArr) ./ pArr;
plot(x, y)
xlabel('Frequency (Hz)')
ylabel('Power (watts)')
title('Frequency vs. Power')
grid on;
axis([-20e9 20e9 -20e9 20e9]);
我希望这段代码能够 运行 顺利地运行,而不会让我的处理器陷入困境。但是,它 运行 并不顺利,并且在我的计算机上用完了两个线程。
bandpass
函数直接过滤提供的输入样本,从而在同一域中生成 'smoothed' 版本。在你的情况下,输入已经是频域表示,所以你会得到一个平滑的频域表示,这不是你想要实现的(这将是平滑时域表示的力量).
相反,由于您已经有了频域表示,您可以简单地 select 正确的频率仓(这基本上就是在频域中应用矩形 window 所做的) .
p0 = sum(pNew((floor(1*n/FS)+1):(floor(60*n/FS)+1)));
p1 = sum(pNew((floor(60*n/FS)+1):(floor(250*n/FS)+1)));
p2 = sum(pNew((floor(250*n/FS)+1):(floor(2e3*n/FS)+1)));
p3 = sum(pNew((floor(2e3*n/FS)+1):(floor(8e3*n/FS)+1)));
p4 = sum(pNew((floor(8e3*n/FS)+1):(floor(20e3*n/FS)+1)));
如果您仍要对 FFT 进行下采样,您可能希望等效地(但更有效地)对较小的 FFT 求和:
% pad to a size that is a multiple of the block size
D = 10;
L = ceil(n/D);
M = L*D;
s2 = zeros(1,M);
s2(1:length(s1)) = s1;
% compute the downsampled FFT by summing smaller FFTs
p = zeros(1,L);
for i = 0:(D-1)
p = p + fft(s2((1+i*L):(1+L+i*L)));
end
问题
我接到一项任务,要在 matlab 中使用 FFT 设计一个音频可视化工具。我已经完成了大部分我想用音频可视化工具做的事情,但我 运行 正在顺利地解决 运行 程序的问题,因为它被我的处理器限制了。我想用我的可视化工具做更多的事情,但我不想再实现任何功能,然后才能很好地掌握 运行 的基础知识。
背景
该程序实际上分为两部分。第一部分是音频播放器和 运行 一个名为 CLICK_ME.m
的 matlab 脚本。第二个是currentfft()
。它从歌曲中获取信息,在我厌倦它之前,它绘制了(频率,功率)。然而,这是最基本的可视化工具,我的教授希望我们创建一个独特的音频可视化工具。
我的想法是绘制其中一个统计五边形,其中五个点中每个点的半径是在五个频段之间过滤的功率。
自然地,我实现了 matlab 的 bandpass()
函数,这实现了它的设计目的。然而,它需要的处理能力比我预期的要多得多,随后它遇到了瓶颈,看起来真的很不稳定。
我尝试通过使用 for 循环每隔一个或每隔三个值来减少必须执行计算的点数。我已经这样做了,直到商数减少了 10 为止,但没有取得明显的成功。如果我删除超过十分之一的点,我会从 matlab 中得到错误,它不能对少于六个样本进行 bandpass()
。
我的教授建议使用 gaussmf()
或 rectangularPulse()
来模拟接近带通滤波器的东西,但我不确定如何实现它们或者它们是否会更快。
function currentfft ( player, Y, FS )
sampleNumber = get( player, 'CurrentSample' );
timerVal = get( player, 'TimerPeriod' );
%Get channel one values for our window around the current sample number
s1 = Y(floor(sampleNumber-((timerVal*FS)/2)):floor(sampleNumber+((timerVal*FS)/2)),1);
n = length(s1);
p = fft(s1); % take the fourier transform
nUniquePts = ceil((n+1)/2);
p = p(1:nUniquePts); % select just the first half since the second half
% is a mirror image of the first
p = abs(p); % take the absolute value, or the magnitude
p = p/n; % scale by the number of points so that
% the magnitude does not depend on the length
% of the signal or on its sampling frequency
p = p.^2; % square it to get the power
% multiply by two
if rem(n, 2) % odd nfft excludes Nyquist point
p(2:end) = p(2:end)*2;
else
p(2:end -1) = p(2:end -1)*2;
end
% reduce the number of points actually being filtered
q = 1;
s = 1;
d = int16( length(p) );
pNew = zeros( [ d/10, 1 ] );
while q < d
pNew(s) = p(q);
q = q + 10;
s = s + 1;
end
% try gaussian or rect functions instead of bandpass?
% radius of each section of the pentagon
p0 = abs( bandpass(pNew, [ 1 60 ], FS) );
p1 = abs( bandpass(pNew, [ 60 250 ], FS) );
p2 = abs( bandpass(pNew, [ 250 2e3 ], FS) );
p3 = abs( bandpass(pNew, [ 2e3 8e3 ], FS) );
p4 = abs( bandpass(pNew, [ 8e3 20e3 ], FS) );
% length( p0 )
% length( p1 )
% length( p2 )
% length( p3 )
% length( p4 )
pArr = [ p0, p1, p2, p3, p4 ];
%freqArray = (0:nUniquePts-1) * (FS / n); % create the frequency array
thetaArr = [ pi/2, 4.5*pi/5, 6.5*pi/5, 8.5*pi/5, 10.5*pi/5 ];
% calculating x and y from the radii
x = cos(thetaArr) ./ pArr;
y = sin(thetaArr) ./ pArr;
plot(x, y)
xlabel('Frequency (Hz)')
ylabel('Power (watts)')
title('Frequency vs. Power')
grid on;
axis([-20e9 20e9 -20e9 20e9]);
我希望这段代码能够 运行 顺利地运行,而不会让我的处理器陷入困境。但是,它 运行 并不顺利,并且在我的计算机上用完了两个线程。
bandpass
函数直接过滤提供的输入样本,从而在同一域中生成 'smoothed' 版本。在你的情况下,输入已经是频域表示,所以你会得到一个平滑的频域表示,这不是你想要实现的(这将是平滑时域表示的力量).
相反,由于您已经有了频域表示,您可以简单地 select 正确的频率仓(这基本上就是在频域中应用矩形 window 所做的) .
p0 = sum(pNew((floor(1*n/FS)+1):(floor(60*n/FS)+1)));
p1 = sum(pNew((floor(60*n/FS)+1):(floor(250*n/FS)+1)));
p2 = sum(pNew((floor(250*n/FS)+1):(floor(2e3*n/FS)+1)));
p3 = sum(pNew((floor(2e3*n/FS)+1):(floor(8e3*n/FS)+1)));
p4 = sum(pNew((floor(8e3*n/FS)+1):(floor(20e3*n/FS)+1)));
如果您仍要对 FFT 进行下采样,您可能希望等效地(但更有效地)对较小的 FFT 求和:
% pad to a size that is a multiple of the block size
D = 10;
L = ceil(n/D);
M = L*D;
s2 = zeros(1,M);
s2(1:length(s1)) = s1;
% compute the downsampled FFT by summing smaller FFTs
p = zeros(1,L);
for i = 0:(D-1)
p = p + fft(s2((1+i*L):(1+L+i*L)));
end