Передискретизируйте простую линейную последовательность в 3/2 исходный уровень 10 Гц. Постройте исходные и передискретизируемые сигналы на той же фигуре.

fs = 10;
t1 = 0:1/fs:1;
x = t1;
y = resample(x,3,2);
t2 = (0:(length(y)-1))*2/(3*fs);

plot(t1,x,'*',t2,y,'o')
xlabel('Time (s)')
ylabel('Signal')
legend('Original','Resampled', ...
    'Location','NorthWest')

При фильтрации resample принимает, что входная последовательность, x, является нулем до и после выборок, которые это дано. Большие отклонения от нуля в конечных точках x могут привести к неожиданным значениям для y.

Покажите эти отклонения путем передискретизации треугольной последовательности и вертикально переключенной версии последовательности с ненулевыми конечными точками.

x = [1:10 9:-1:1;
    10:-1:1 2:10]';
y = resample(x,3,2);

subplot(2,1,1)
plot(1:19,x(:,1),'*',(0:28)*2/3 + 1,y(:,1),'o')
title('Edge Effects Not Noticeable')
legend('Original','Resampled', ...
    'Location','South')

subplot(2,1,2)
plot(1:19,x(:,2),'*',(0:28)*2/3 + 1,y(:,2),'o')
title('Edge Effects Noticeable')
legend('Original','Resampled', ...
    'Location','North')

Создайте синусоидальный сигнал. Задайте частоту дискретизации, таким образом, что 16 выборок соответствуют точно одному периоду сигнала. Чертите диаграмму стебель-листья сигнала. Наложите ступенчатый график для визуализации выбирать-и-содержать.

fs = 16;
t = 0:1/fs:1-1/fs;

x = 0.75*sin(2*pi*t);

stem(t,x)
hold on
stairs(t,x)
hold off

Используйте resample, чтобы сверхдискретизировать сигнал фактором четыре. Используйте настройки по умолчанию. Постройте результат вместе с исходным сигналом.

ups = 4;
dns = 1;

fu = fs*ups;
tu = 0:1/fu:1-1/fu;

y = resample(x,ups,dns);

stem(tu,y)
hold on
stairs(t,x)
hold off
legend('Resampled','Original')

Повторите вычисление. Задайте n = 1 так, чтобы фильтр сглаживания имел порядок $$2\times 1\times 4=8$$. Задайте параметр формы $$\beta =0$$ для окна Kaiser. Выведите фильтр, а также передискретизируемый сигнал.

n = 1;
beta = 0;

[y,b] = resample(x,ups,dns,n,beta);

fo = filtord(b)

fo = 8

stem(tu,y)
hold on
stairs(t,x,'--')
hold off
legend('n = 1, \beta = 0')

Передискретизируемый сигнал показывает эффекты искажения, которые следуют из относительно широкого mainlobe и низкого затухания бокового лепестка окна.

Увеличьте n до 5 и отпуск $$\beta =0$$. Проверьте, что фильтр имеет порядок 40. Постройте передискретизируемый сигнал.

n = 5;

[y,b] = resample(x,ups,dns,n,beta);

fo = filtord(b)

fo = 40

stem(tu,y)
hold on
stairs(t,x,'--')
hold off
legend('n = 5, \beta = 0')

Более длинное окно имеет более узкий mainlobe и ослабляет эффекты искажения лучше. Это также ослабляет сигнал.

Оставьте порядок фильтра в $$2\times 5\times 4=\mathrm{}40$$ и увеличьте параметр формы до $$\beta =\mathrm{}20$$.

beta = 20;

y = resample(x,ups,dns,n,beta);

stem(tu,y)
hold on
stairs(t,x,'--')
hold off
legend('n = 5, \beta = 20')

Высокое затухание бокового лепестка приводит к хорошей передискретизации.

Уменьшите порядок фильтра назад к $$2\times 1\times 4=8$$ и отпуск $$\beta =\mathrm{}20$$.

n = 1;

[y,b] = resample(x,ups,dns,n,beta);

stem(tu,y)
hold on
stairs(t,x,'--')
hold off
legend('n = 1, \beta = 20')

Шире mainlobe генерирует значительные артефакты после передискретизации.

Сгенерируйте 60 выборок синусоиды и передискретизируйте ее в 3/2 исходный уровень. Отобразите исходные и передискретизируемые сигналы.

tx = 0:6:360-3;
x = sin(2*pi*tx/120);

ty = 0:4:360-2;
[y,by] = resample(x,3,2);

plot(tx,x,'+-',ty,y,'o:')
legend('original','resampled')

Постройте частотную характеристику фильтра сглаживания.

freqz(by)

Передискретизируйте сигнал в 2/3 исходный уровень. Отобразите исходный сигнал и его передискретизацию.

tz = 0:9:360-9;
[z,bz] = resample(x,2,3);

plot(tx,x,'+-',tz,z,'o:')
legend('original','resampled')

Постройте импульсный ответ нового фильтра lowpass.

impz(bz)

Используйте данные, зарегистрированные Галилео Галилеем в 1 610, чтобы определить орбитальный период Каллисто, наиболее удаленный из четырех самых больших спутников Юпитера.

Галилео наблюдал движение спутников в течение шести недель, начиная 15 января. Наблюдения имеют несколько разрывов, потому что Юпитер не был видим облачными ночами. Сгенерируйте массив datetime времен наблюдения.

t = [0 2 3 7 8 9 10 11 12 17 18 19 20 24 25 26 27 28 29 31 32 33 35 37 ...
    41 42 43 44 45]'+1;

yg = [10.5 11.5 10.5 -5.5 -10.0 -12.0 -11.5 -12.0 -7.5 8.5 12.5 12.5 ...
    10.5 -6.0 -11.5 -12.5 -12.5 -10.5 -6.5 2.0 8.5 10.5 13.5 10.5 -8.5 ...
    -10.5 -10.5 -10.0 -8.0]';

obsv = datetime(1610,1,15+t);

Передискретизируйте данные на обычную сетку с помощью частоты дискретизации одного наблюдения в день. Используйте умеренный фактор повышающей дискретизации 3, чтобы не сверхсоответствовать.

fs = 1;

[y,ty] = resample(yg,t,fs,3,1);

Отобразите на графике данные и передискретизируемый сигнал.

plot(t,yg,'o',ty,y,'.-')
xlabel('Day')

Повторите процедуру с помощью интерполяции сплайна и отобразив даты наблюдения. Выразите частоту дискретизации в обратные дни.

fs = 1/86400;

[ys,tys] = resample(yg,obsv,fs,3,1,'spline');

plot(t,yg,'o')
hold on
plot(ys,'.-')
hold off

ax = gca;
ax.XTick = t(1:9:end);
ax.XTickLabel = char(obsv(1:9:end));

Вычислите оценку спектра мощности периодограммы однородно расположенных с интервалами, линейно интерполированных данных. Выберите длину DFT 1 024. Сигнал достигает максимума при инверсии орбитального периода.

[pxx,f] = periodogram(ys,[],1024,1,'power');
[pk,i0] = max(pxx);

f0 = f(i0);
T0 = 1/f0

T0 = 16.7869

plot(f,pxx,f0,pk,'o')
xlabel('Frequency (day^{-1})')