Вычисление ДПФ произвольного сигнала

Question

В рамках курса по обработке сигналов в университете нас попросили написать алгоритм в Matlab для расчета одностороннего спектра нашего сигнала с использованием DFT, без использования функции fft(), встроенной в matlab.это не оцениваемая часть курса, я просто заинтересован в том, чтобы сделать это «правильным» для себя.В настоящее время я использую версию Matlab 2018b, если кто-нибудь сочтет это полезным.

Я построил сигнал синусоиды 1 кГц и 2 кГц, сдвинутый по фазе на 135 градусов (2 * пи / 3 рад).

затем, используя уравнения в 9.1 из Обработка сигнала с дискретным временем (Аллан В. Оппенгейм) и формулу Эйлера для упрощения показателя степени, я создаю этот код:

%%DFT(currently buggy)
n=0;m=0;
for m=1:DFT_N-1 %DFT_Fmin;DFT_Fmax; %scrolls through DFT m values (K in text.)
    for n=1:DFT_N-1;%;(DFT_N-1);%<<redundant code? from Oppenheim eqn. 9.1 % eulers identity, K=m and n=n
        X(m)=x(n)*(cos((2*pi*n*m)/DFT_N)-j*sin((2*pi*n*m)/DFT_N));
        n=n+1;
    end
    %m=m+1; %redundant code?
end

Это принимает х в качестве входа, в этом случае сигнал, упомянутый ранее, а также разрешение преобразования, как представлено DFT_N, который был инициализирован до 100. Выход этой функции, X, должен быть чем-то вчастотная область, но построение графика X дает круговой график, немного больший, чем единичный круг, и с пропуском по левому краю.

Я изо всех сил пытаюсь понять, как я должен преобразовать это в stem()графики, представленные встроенным алгоритмом DFT.

Большое спасибо, J.

Answer 1

Если вас интересует только односторонний спектр, то вы можете просто рассчитать X (m) for m=1:DFT_N/2, поскольку X сопряженно симметричен относительно m = DFT_N / 2, т. Е. X(DFT_N/2+m) = X(DFT_N/2-m)', из-за exp(-j*(pi*n+2*pi/DFT_N*m)) = exp(-j*(pi*n-2*pi/DFT_N*m))'.

В качестве примечания, для данной m эта программа вычисляет внутреннее произведение между массивом x и другим массивом комплексных экспонент, т. Е. exp(-j*2*pi/DFT_N*m*n), для n = 0,1, ..., N-1.Синтаксис MATLAB очень удобен для таких вычислений, и вы можете избежать этого внутреннего цикла с помощью следующей команды

exp(-j*2*pi/DFT_N*m*(0:DFT_N-1)) * x

, где x - вектор столбца.Точно так же вы можете избежать первого цикла, расширяя ваш комплексный экспоненциальный вектор по строкам для каждого m, т. Е. Строить матрицу exp(-j*2*pi/DFT_N*(0:DFT_N-1)'*(0:DFT_N-1)).Тогда ваш ДПФ просто

X = exp(-j*2*pi/DFT_N*(0:DFT_N-1)'*(0:DFT_N-1)) * x

Для одностороннего спектра, вместо этого используйте

X = exp(-j*2*pi/DFT_N*(0:floor((DFT_N-1)/2))'*(0:DFT_N-1)) * x

Answer 2

Это ваша ошибка:

заменить X(m)=x(n)*(cos.. на X(m)=X(m)+x(n)*(cos..

Для заданного m он не интегрируется по переменной n, а перезаписывает X (m) толькоПоследнее вычисление для n = DFT_N-1.

Обратите внимание, что при интегрировании по n=1:DFT_N-1 пропускается одна гармоника, т. е. первая гармоника exp (-j * 2 * pi).Замените n=1:DFT_N-1 на n=1:DFT_N, чтобы включить это.Я бы также заменил m=1:DFT_N-1 на m=1:DFT_N для составления графика проблем.

Также заменим любой 2*pi*n*m на 2*pi*(n-1)*(m-1), чтобы получить правильную фазу, поскольку первая запись X должна соответствовать нулевой частоте,получая sum_n x (n) * (cos (0) + j sin (0)) = sum_n x (n).Если ваш сигнал x является действительным значением, тогда компонент нулевой частоты X (1) должен быть действительным значением, угол (X (1)) = 0.

Последнее замечание, не забудьте сдвинуть ноль-частотная составляющая к центру спектра для лучшей видимости, X = circshift(X,floor(size(X)/2));