#matlab #signal-processing #fft
#matlab #обработка сигналов #БПФ
Вопрос:
Проблема
Мне было дано задание разработать аудиовизуализатор с использованием FFT в matlab. Я прошел большую часть того, что я хочу сделать с аудиовизуализатором, но у меня возникают проблемы с плавным запуском программы, потому что она регулируется моим процессором. Я хочу сделать больше с моим визуализатором, но я не хочу реализовывать больше функций, прежде чем я смогу заставить основы работать хорошо.
Фон
Программа на самом деле состоит из двух частей. Первая часть — это аудиоплеер, который запускает вызываемый скрипт matlab CLICK_ME.m . Второй currentfft() . Он принимает информацию из песни и, прежде чем я ее загрузил, отображает (частоту, мощность). Однако это самый простой визуализатор, и мой профессор хочет, чтобы мы создали уникальный аудиовизуализатор.
Моя идея заключалась в том, чтобы отобразить один из этих пятиугольников, где радиус каждой из пяти точек равен мощности, фильтруемой между пятью полосами частот.
Естественно, я реализовал bandpass() функцию matlab, и она выполняет то, для чего она предназначена. Однако для этого требуется намного больше вычислительной мощности, чем я мог ожидать, и впоследствии он становится узким местом и выглядит действительно прерывистым.
Я попытался уменьшить количество точек, с которыми он должен выполнять вычисления, используя цикл for для получения каждого второго или каждого третьего значения. Я сделал это с уменьшением на коэффициент десять без заметного успеха. Если я удаляю более каждой десятой точки, я получаю ошибки от 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]);
Я бы хотел, чтобы этот код работал плавно, не пропуская мой процессор через отжим. Однако он не работает гладко и использует максимум двух потоков на моем компьютере.
Комментарии:
1. Я не понимаю, зачем вообще нужны полосовые фильтры. Если у вас уже есть частотный спектр, вы можете просто рассмотреть отдельные интересующие вас полосы вашего спектра. Таким образом, это все равно что применить идеальный полосовой фильтр к вашему сигналу. Просто возьмите первые 60 бинов вашего спектра для первой полосы, следующие 60 бинов для второй полосы и так далее.
2. Просто для уточнения. Этот подход сводился бы к умножению вашей частотной характеристики (вывода БПФ) на сдвинутое прямоугольное окно для получения отдельной полосы
3. К сожалению, мы еще не перешли к сдвинутым прямоугольным окнам. Я знаю, что существует функция с именем rectwin(), которая принимает один аргумент для длины окна, но я не уверен, как его можно сдвинуть? Я просматривал документацию по обработке звука и обработке сигналов, чтобы найти что-нибудь полезное, но безуспешно.
Ответ №1:
bandpass Функция напрямую фильтрует предоставленные входные выборки, что приводит к «сглаженной» версии в том же домене. В вашем случае этот ввод уже является представлением в частотной области, поэтому вы получите сглаженное представление в частотной области, чего вы не пытаетесь достичь (вместо этого это будет сила сглаженного представления во временной области).
Вместо этого, поскольку у вас уже есть представление в частотной области, вы можете просто выбрать правильные частотные ячейки (что, по сути, и делает применение прямоугольного окна в частотной области).
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)));
Если вы все еще собираетесь понижать дискретизацию вашего БПФ, вы можете вместо этого захотеть эквивалентно (но более эффективно) суммировать меньшие БПФ:
% 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