Я пытаюсь написать простую функцию в Matlab для расчета и построения спектра мощности шума (NPS). Прежде всего, я хотел проверить, подходит ли алгоритм, который я получил от моего учителя, и все такое. Вот оно (сделано в mathcad)
Поэтому я попытался скопировать и вставить его в скрипт Matlab и в итоге получил этот код:
clear all;
clc;
N=1000;
O=1024;
mn=zeros(N,O);
n0=1500;
s=sqrt(n0);
W=zeros(N,O/2);
W1=zeros(N,O);
for k=1:N
for l=1:O
mn(k,l)=n0+round(sin(randn)*s);
end
end
for k=1:N
for l=1:O
mn(k,l)=mn(k,l)-n0;
end
end
for k=1:N
W1(k,:)=fft(mn(k,:));
end
for k=1:N
for l=1:O/2
W(k,l)=W1(k,l);
end
end
NPS1=(abs(W)).^2;
NPS2=sum(NPS1);
NPS=(1/N)*NPS2;
plot(NPS);
Я не использую распределение Пуассона, и я переключил индексы строка-столбец, но это не имеет значения (верно?). Проблема в том, что значения на моем графике в 400 раз больше, чем ожидалось.
Вот как это должно выглядеть:
Я пытался найти, что я сделал не так, но после некоторого времени тестирования некоторых изменений я вернулся в исходное положение ... Единственное, что меня беспокоит, это то, что, возможно, функция Matlab fft работает иначе, чем та, которая использовалась в Mathcad (не могу сказать, что полностью понимаю). Любая добрая душа может сказать мне, если дело в функции fft? Или я просто слепой болван, который не видит совершенную им глупую ошибку? Приветствую и извините за мой плохой английский.
[РЕДАКТИРОВАТЬ]
Через некоторое время мой учитель попросил меня проверить, работает ли этот метод с коррелированным (своего рода) шумом, поскольку он работает (снова) в mathcad. После корреляции его NPS должен «падать» на более высоких частотах. Проблема в том, что это не так. Код, который я использую для проверки, выглядит так:
clear all;
clc;
N=1000;
mn = poissrnd(N, N, N);
dataw=zeros(N);
for k=1:N ## loop used for my teacher's correlation method
for l=1:N
if l>1 && l<N
dataw(k,l)=dataw(k,l)+mn(k,l)*0.5+mn(k,l-1)*0.25+mn(k,l+1)*0.25;
elseif l==1
dataw(k,l)=dataw(k,l)+mn(k,l)*0.75+mn(k,l+1)*0.25;
else
dataw(k,l)=dataw(k,l)+mn(k,l)*0.75+mn(k,l-1)*0.25;
end
end
end
dataw = dataw - mean(dataw(:));
W1 = (1/sqrt(N))*fft(dataw, [], 1);
NPS1=(abs(W1)).^2;
NPS2=sum(NPS1);
NPS=(1/N)*NPS2;
plot(NPS);
Единственные изменения, которые я внес в код, исправленный Rayryeng, - это квадрат матрицы шума (1000x1000), а также среднее значение 1000 и использование всего преобразованного вектора W1 вместо его «половины». Я знаю, что это сработало для моего учителя, но это не для меня ... Есть ли что-то еще в Matlab fft, которое я упустил, или это «метод корреляции», который я использовал?
Добавление того, как это выглядит в Mathcad после нескольких быстрых изменений (несколько незначительных отличий, но в целом это показывает эффект, который я должен получить). Это прервало начало сканирования, но оно точно такое же, как я указал в начале этого поста.
[РЕДАКТИРОВАТЬ2]
Nvm, это была проблема размера в функции fft. После изменения на fft (dataw, [], 2) он выглядит лучше.
