Изменение фазы сигнала в частотной области (MatLab)

Question

Я отправил этот вопрос на dsp.stackexchange, и был проинформирован о том, что это более актуально для StackOverflow, как это в первую очередь программирования вопрос:

Я пытаюсь написать код, который позволяет мне изменять фазу сигнала в частотной области. Однако мой вывод не совсем корректен, поэтому что-то должно быть неправильным. Для простого примера предположим, что мы имеем функцию y = sin (2 * pi * t) и хотим реализовать фазовый сдвиг -pi/2. Мой код выглядит следующим образом:

clear all 
close all 

N = 64; %number of samples 
fs = 10; %sampling frequency 
ts = 1/fs; %sample interval 
tmax = (N-1)*ts; 
t = 0:ts:tmax; 
y = sin(2*pi*t); 

figure 
plot(t,y) 

% We do the FT 
f = -fs/2:fs/(N-1):fs/2; 
Y = fftshift(fft(y)); 

% Magnitude spectrum 
figure 
plot(f,abs(Y)); 

phase = angle(Y); 

% Phase spectrum 
figure 
plot(f,phase) 

Y = ifftshift(Y) 

% Attempt at phase shift 
Y = Y.*exp(-i*2*pi*f*pi/2); 

% Inverse FT 
u = ifft(Y); 

figure 
plot(t,real(u)) 

Все участки исправны для конечного участка, который выглядит следующим образом, за исключением:

Это выглядит почти правильно, но не совсем. Если кто-нибудь может дать мне несколько указаний относительно того, как мой код может быть исправлен, чтобы исправить это, я был бы очень признателен! У меня такое чувство, что моя ошибка связана с линией Y = Y.*exp(-i*2*pi*f*pi/2);, но я не уверен, как ее исправить.

Answer 1

Я не могу вдаваться в подробности Фурье анализа (потому что я не знаю их), но я могу предложить рабочее решение с некоторыми намеками:

Прежде всего, Вы должны выразить вашу волну мнимые члены, а именно:

y = exp(1i*2*pi*t); 

И что еще более важно, Вы должны действительно перейти только фазу, не мешая со всем спектром:

% Attempt at phase shift 
Y = abs(Y).*exp(1i*angle(Y)-1i*pi/4); % -pi/4 shift 

Вы шо uld отмечают, что сдвиг больше не связан с частотой, что, на мой взгляд, имеет смысл. Наконец Вы можете начертить результаты:

figure 
plot(t,real(u),'k') 
hold on 
plot(t,real(y),'r') 

real(y) фактически является функцией косинуса, и Вы начали с синусоидальной, но, надеюсь, вы получите идею. Для пи/4 сдвига я получил что-то вроде этого (начался с красным, отделанным чернымом):

Answer 2

Вы сделал 3 главных ошибки в дизайне коды.

1. Входной вектор БПФ интерпретируется как период сигнала, который повторяется бесконечно. Это означает, что ваш входной вектор должен содержать целое число полных периодов вашего синусоидального сигнала . У вас есть входной вектор из 64 выборок и частота выборки 10. Это приводит к 6,4 периодам вашей синусоидальной волны, что приводит к утечке. Если вы проверите частотный спектр после , выполняя БПФ, вы увидите, что нет двух чистых частотных линий, но много частотных компонентов вокруг двух мест.
2. После исправления вашего входного вектора должны быть только две одиночные частоты со значениями , которые не близки к нулям. 62 будут состоять из численного шума очень , близкого к нулю. Вычисление фазы этих значений приводит к данным мусора.
3. Фазовый сдвиг pi/2 во временной области эквивалентен сдвигом во временной области на N/4, если N равно количеству входных выборок.

Я изменил ваш код. Вы найдете его ниже. С переменной M вы можете изменить количество периодов синусоидальной волны в вашем входном векторе. В примере я установил M = 3.

clear all; 
close all; 

T = 1; %length of sampling sequence in s 
N = 64; %number of samples 
M = 3; % number of periods per sequence 
ts = T/N; %sample interval 
fs = 1/ts %sampling frequency 
tmax = (N-1)*ts; 
t = 0:ts:tmax; 
y = sin(2*pi*M*t); 

fig01 = figure; 
plot(t,y); 
grid on; 

%% We do the FT 
Y = fft(y); 

%% We create a frequency vector in natural order 
% -fs/2, ..., 0, ... +fs/2 
f =fftshift((0:(fs-1)) - fs/2); 

%% Show Magnitude spectrum 
% There shold be only two lines at -M and +M 
figure; 
plot(f,abs(Y),'o'); 
grid on; 

%% Attempt at phase shift 
% y/t) -> Y(w) 
% y(t-t0) -> Y(w) * exp(-i*w*t0) 
% Phase shift of pi/2 in frequncy domain is equavalent to as time shift 
% of T/4 in time domain 

Y = Y.*exp(-i*2*pi*f*T/4); 

% Inverse FT 
u = ifft(Y); 

figure 
hold on; 
plot(t,real(u),'b-'); 
plot(t,real(y),'r-'); 
hold off; 
grid; 

Входной сигнал с трех периодов синусоидального сигнала

спектра входного сигнала. Частотные линии на -3 и +3

Входной сигнал (синий) и сдвинуты по фазе сигнал (красный)