Здесь я хочу расписать некоторые аспекты реализации дискретного преобразования Фурье (DFT)
Сначала я крайне рекомендую почитать вот этот блогопост. Очень толково и очень просто расписано про ключевые моменты DFT.
Несмотря на степень разжевывания в той статье, у меня все равно оставились некоторые вопросы, которые я попытался раскопать и порешать.
Идем по-порядку, сначала нужно понимать что такое корреляция – это по сути сумма продуктов умножения двух функций друг на друга. Может быть как положительной так и отрицательной
\sum_{i=0}^{N}x(i)y(i)Говорим про обычную корреляцию, без всяких data science штук 🙂
Теперь примерно вырисовывается базовая идея вложенная в дискретное преобразование Фурье – мы берем нашу функцию, умножаем ее на известную нам простую функцию с одной гармоникой и смотрим на каких именно гармониках (или бинах) у нас корреляция больше и строим таким образом спектр.
Как обычно, самый простой и понятный способ это пройти через какой-нибудь пример, предположим что у нас есть синус и косинус на 7Гц.
Немного кода в Octave:
clf
clear
#input data variables
step = 1/256;
tstop = 2;
f1 = 7;
t = 0:step:tstop-step;
for i = 1:length(t)
rand_coeff1(i) = (rand()-0.5:0.5)/3;
rand_coeff2(i) = (rand()-0.5:0.5)/3;
end
y = sin(2*pi*f1*t);
ytry=e.^(j*2*pi*7*t);
y_product1 = y.*ytry;
plot(t,y)
hold on
grid on
correlation1 = sum(imag(y_product1));
correlation2 = sum(real(y_product1));
y = sin(2*pi*f1*t+pi/2);
plot(t,y)
ytry=e.^(j*2*pi*7*t);
y_product1 = y.*ytry;
correlation3 = sum(imag(y_product1));
correlation4 = sum(real(y_product1));
Я в коде не использовал магнитуды, иначе мы не сможем отделить синус от косинуса.
| correlaton1 | 256 |
| correlation2 | ~10^(-14) |
| correlation3 | ~10^(-13) |
| correlation4 | 256 |
Отсюда видим что корреляция1 которая берет мнимую часть равна 256, а корреляция 4 берет реальную часть, а значит 1 это синус, а 4 это косинус.
Теперь можно вернуться к нашему примеру, немного его усложнить – теперь у нас пускай существует два сигнала – на 7 и на 8Гц.
clf
clear
#input data variables
step = 1/256;
tstop = 2;
f1 = 7;
t = 0:step:tstop-step;
N = 16;
y = sin(2*pi*f1*t)+sin(2*pi*(f1+1)*t);
#ytry=e.^(j*2*pi*1*t);
#y_product1 = y.*ytry;
#plot(t,y)
##hold on
#grid on
for i=1:1:N
ytry=e.^(j*2*pi*i*t);
y_product = y.*ytry;
correlation_real(i) = sum(real(y_product));
correlation_imag(i) = sum(imag(y_product));
correlation_abs(i) = correlation_real(i) + correlation_imag(i);
plot(i,correlation_abs(i),'marker','o');
hold on;
grid on;
endfor
Проверяем чего вышло:
Здесь отчетливо видны два пика на 7 и на 8 бине, но сразу вылезает пара вопросов: Почему мы получаем огромное значение по y? Как правильно опредить N? Какие лимиты должны быть и как вообще нормально сделать график 🙂
Переходим к более реалистичной функции, и сейчас используем это все вместе с цифрами, будет за что нам зацепиться.
y = 2 + 5cos(2 \pi2 t)+3sin(2\pi3t)+2cos(2\pi4t)
Итак синяя кривая описана функцией, приведенной выше, в то время как красная это сэмплированные данные с частотой Fs = 8Гц, и я засемплил 8 точек, значит наше разрешение 1Гц.
Попытаемся сначала ответить на вопрос почему у нас по Y такие большие числа. Очевидно что это как то соотносится с количеством точек, которые были использованы для расчета корреляции. Первое желание – просто разделить все результаты на N, посмотрим что из этого получится:
clf
clear
#input data variables
step = 1/1024;
tstop = 1;
t = 0:step:tstop;
y = 2 + 5*cos(2*pi*t*2) + 3*sin(2*pi*t*3) + 2*cos(2*pi*t*4);
#plot(t,y,'b', 'linewidth',4);
grid on;
hold on;
N = 8;
Fs = 8;
ts = 1/Fs;
f_spectrum = 0:Fs/N:(Fs*(N-1))/N;
tsample = 0:ts:(N-1)*ts;
y_sampled = 2 + 5*cos(2*pi*tsample*2)+3*sin(2*pi*tsample*3) + 2*cos(2*pi*tsample*4);
#plot(tsample,y_sampled,'r','marker','o','linewidth',1,'markersize',20);
h=legend('original ','sampled ');
set(h,'fontsize',14,'fontname','FreeSans','fontweight','normal');
set(gca, 'linewidth', 3, 'fontsize', 18,'fontname', 'FreeSans') #modify line and fontsize of x and y axes
#xlabel('t, s');
#ylabel('Amplitude, V');
for i=0:1:(N-1)
ytry=e.^(-j*2*pi*i*Fs/N*tsample);
y_product = (y_sampled.*ytry)/N;
correlation_real(i+1) = sum(real(y_product));
correlation_imag(i+1) = sum(imag(y_product));
correlation_abs(i+1) = sqrt(correlation_real(i+1)^2 + correlation_imag(i+1)^2);
#plot(i,correlation_abs(i+1),'marker','o');
#hold on;
#grid on;
endfor
xlabel('F, Hz');
ylabel('Amplitude, V');
bar(f_spectrum,correlation_abs,0.1)
Так, пару вопросов в принципе отвечено здесь, но есть пару НО. Первое – величина гармоники по DC и на Fs/2 (0 и 4Hz) соответствует правильным, а вот все что внутри промежутка равно половине. Также мы видим сигналы после 4Гц, что есть алиасинг, мы вошли во вторую зону Найквиста и получили “складывание” спекта (только во второй половине все сигналы комплексно сопряженные). Например если у нас 6Гц сигнал наблюдается, это просто 8-2 отражение. Именно это и есть причина деградации спектра в первой зоне Найквиста, в то время как точки в 0 и 4 совпадают и со “сложенным” спектром. Отсюда можно сделать вывод -все точки внутри DC – Fs/2 надо поделить на N/2 а не на N. Также для dft построения спектра, нас в принципе не особо интересуют точки после Fs/2.
clf
clear
#input data variables
step = 1/1024;
tstop = 1;
t = 0:step:tstop;
y = 2 + 5*cos(2*pi*t*2) + 3*sin(2*pi*t*3) + 2*cos(2*pi*t*4);
#plot(t,y,'b', 'linewidth',4);
grid on;
hold on;
N = 8;
Fs = 8;
ts = 1/Fs;
f_spectrum = 0:Fs/N:(Fs*(N-1))/N;
tsample = 0:ts:(N-1)*ts;
y_sampled = 2 + 5*cos(2*pi*tsample*2)+3*sin(2*pi*tsample*3) + 2*cos(2*pi*tsample*4);
#plot(tsample,y_sampled,'r','marker','o','linewidth',1,'markersize',20);
h=legend('original ','sampled ');
set(h,'fontsize',14,'fontname','FreeSans','fontweight','normal');
set(gca, 'linewidth', 3, 'fontsize', 18,'fontname', 'FreeSans') #modify line and fontsize of x and y axes
#xlabel('t, s');
#ylabel('Amplitude, V');
for i=0:1:(N-1)
ytry=e.^(-j*2*pi*i*Fs/N*tsample);
if ((i==0)||(i==N/2))
y_product = (y_sampled.*ytry)/N;
else
y_product = (y_sampled.*ytry)/(0.5*N);
endif
correlation_real(i+1) = sum(real(y_product));
correlation_imag(i+1) = sum(imag(y_product));
correlation_abs(i+1) = sqrt(correlation_real(i+1)^2 + correlation_imag(i+1)^2);
#plot(i,correlation_abs(i+1),'marker','o');
#hold on;
#grid on;
endfor
subplot(3,1,1);
title("mag");
bar(f_spectrum,correlation_abs,0.1)
subplot(3,1,2);
title("real")
bar(f_spectrum,correlation_real,0.1)
subplot(3,1,3)
title("imag")
bar(f_spectrum,correlation_imag,0.1)
С таким кодом у нас в принципе получилось получить правильные амплитуды и показать синусы и косинусы, утверждение про комплексную сопряженность хорошо видно по нижнему графику для синуса.
Это все замечательно, но если посмотреть на формулы dft:
F[n] = \sum_{k=0}^{N-1}f[k]e^{-j2\pi kn/N} (n=0:N-1)То можно заметить что она отличается от реализации которую я применил в моем коде. Octave или matlab позволяет перемножать кривые друг с другом, в то время как не все языки программирования имеют эту роскошь (наверное). Вместо этого мы чекаем каждый бин (n) вычисляя корреляционную сумму сигнала с сигналом kn/N. Преобразуем наш код, чтобы привести его в соответствие с оригинальной формулой:
clf;
clear;
grid on;
hold on;
N = 8;
Fs = 8;
ts = 1/Fs;
f_spectrum = 0:Fs/N:(Fs*(N-1))/N;
tsample = 0:ts:(N-1)*ts;
y_sampled = 2 + 5*cos(2*pi*tsample*2)+3*sin(2*pi*tsample*3) + 2*cos(2*pi*tsample*4);
#plot(tsample,y_sampled,'r','marker','o','linewidth',1,'markersize',20);
h=legend('original ','sampled ');
set(h,'fontsize',14,'fontname','FreeSans','fontweight','normal');
set(gca, 'linewidth', 3, 'fontsize', 18,'fontname', 'FreeSans') #modify line and fontsize of x and y axes
#xlabel('t, s');
#ylabel('Amplitude, V');
for i=0:1:(N-1)
sum_harm = 0;
for k=0:1:(N-1)
sum_harm = sum_harm + y_sampled(k+1)*e^(-j*2*pi*i*k/N);
endfor
f(i+1) = sqrt(real(sum_harm)^2+imag(sum_harm)^2);
endfor
#normalization procedure
for i=1:1:N
if ((i==1)||(i==(N/2 + 1)))
f_norm(i) = f(i)/N;
else
f_norm(i) = f(i)*2/N;
endif
endfor
bar(f_spectrum,f_norm,0.1)
xlabel('f, Hz');
ylabel('Amplitude, V');Результат все еще не поменялся:
В принципе с DFT все более менее ясно, осталось пару вещей которые надо взять на заметку:
- Aliasing – это мы уже видели, у нас есть алиаисинг приводит к появлению гармоник после Fs/2
- Spectral leakage – а вот с этим мы еще не сталкивались, в нашем примере все частоты идеально делятся на частоту семплирования и у нас вмещается целочисленное число периодов на нашем промежутке
Например, пускай у нас есть 3.3 Гц частота с амплитудой 3 и частота семплирования все еще 8Гц. Также давайте увеличим выборку N до 64, таким образом разрешение будет 0.125Гц и наша частота не попадает ровно в наше разрешение.
Куча непонятных амплитуд размазана по всему спектру, по всем бинам.
Чтобы уменьшить влияние этого эффекта существует штука под названием окно 🙂 С помощью окна на оригинальный сигнал (до процесса fft) накладывается определенная функция (transfer function) которая делает эффект размытия спектра менее значимым и позволяет выжать больше полезной информации из спектра.
Одно из самых популярных окон это Hanning window:
w(n) = 0.5(1-cos(2\pi \frac{n}{N}))Если ее заплотить то можно увидеть это:
То есть окно делает ни что иное как “сглаживает” края семплированных входных данных ближе к нулю, давайте попробуем реализовать его:
N = 64;
Fs = 8;
ts = 1/Fs;
#f_spectrum = 0:Fs/N:(Fs*(N-1))/N;
f_spectrum = 0:Fs/N:Fs/2;
tsample = 0:ts:(N-1)*ts;
nsample = 0:1:N;
y_sampled = 2 + 3*cos(2*pi*tsample*3.3);
for i=1:1:(N)
ywindowed(i) = 2+ 0.5*(1-cos(2*pi*i/N))*(y_sampled(i)-2);
endfor
plot(tsample,y_sampled,'b','linewidth',2,'markersize',10);
hold on;
plot(tsample,ywindowed,'r','marker','o','linewidth',2,'markersize',10);Я здесь добавил немного пост процессинга сверху чтобы убрать dc компонент а потом добавить, чтобы получить более красивый график:
Если мы теперь возьмем и сделаем dft на данных после того как мы добавили окно:
Конечно все еще не идеальный спектр, но его мы и не получим с ограниченным количеством точек то. Зато теперь нету кучи мусора размытого вообще по всем бинам, вместо того все сконцентрировано вокруг нашей 3.3Гц частоты.
Вообще в октаве или матлабе есть встроенные функции окон, не обязательно свое писать:
ywindowed = y_sampled.*hanning(N)';OK, вроде теперь более мене понятно все. Осталось только упомянуть одну важную штуку – процессорное время DFT. Дело в том что алгоритм дфт требуют просто кучу перемножений – (N*N), хотите больше разрешение получаете больше перемножений, например для 65536 точек в результате вылезет 4294967296 операций, что неплохо так загрузит процессорное время. Чтобы этого избежать и изобрели Fast Fourier Transform (FFT) в котором используется не тупое перемножение, а те коэффициенты которые повторяются просто используются несколько раз. Звучит просто, на деле все намного сложнее и там целая отдельная тема в которую надо копать и вкуривать.
И последнее, давайте просто попробуем использовать встроенную fft функцию вместо того чтобы клепать свои:
clf
clear
#input data variables
tstop = 1;
hold on;
N = 8;
Fs = 8;
ts = 1/Fs;
f_spectrum = 0:Fs/N:(Fs*(N-1))/N;
#define window there
#window = blackman(N);
window = 1;
tsample = 0:ts:(N-1)*ts;
y_sampled = 2 + 5*cos(2*pi*tsample*2)+3*sin(2*pi*tsample*3) + 2*cos(2*pi*tsample*4);
yFFT = y_sampled.*window';
set(gca, 'linewidth', 3, 'fontsize', 18,'fontname', 'FreeSans') #modify line and fontsize of x and y axes
#plot(tsample,y_sampled,'b','linewidth',2,'markersize',20);
grid on;
spectrum_complex = fft(yFFT, N)/(0.5*N);
spectrum_complex(1) = spectrum_complex(1) /2;
spectrum_complex(N/2+1) = spectrum_complex(N/2+1) /2;
spectrum_mag = abs(spectrum_complex);
bar(f_spectrum(1:(N/2+1)), spectrum_mag(1:(N/2+1)),0.2);
Вроде бы похоже на те результаты которые выходили и раньше.