2014-09-21 42 views
0

我张贴在dsp.stackexchange这个问题,并得知它是计算器更加相关,因为它主要是一个编程问题:变化阶段(matlab程序)

我试图写一个代码,它允许我改变频域信号的相位。但是,我的输出不完全正确,所以一定是错误的。举一个简单的例子,假设我们有函数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)) 

所有地块看,除了最后的情节看起来如下OK:

Plot of signal with phase change

这看起来几乎是正确的,但不完全是。如果任何人都可以给我一些关于如何修正我的代码以解决此问题的指示,我将不胜感激!我有一种感觉,我的错误与Y = Y.*exp(-i*2*pi*f*pi/2);行有关,但我不知道如何解决它。

回答

4

我真的不能进入傅立叶分析细节(因为我真的不知道它们),但我可以提供一些提示的工作方案:

首先,你应表示在波虚数项,即:

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)实际上是一个余弦函数,你开始与正弦,但希望你的想法。 对于PI/4移位我有这样的事情(开始用红色,搭配黑色的完成):

this is image description here, how are You?

+0

感谢!这真的很有帮助! – Kristian 2014-09-21 21:36:06

3

你在你的代码设计制作3次重大失误。

  1. FFT的输入矢量被解释为信号的周期,并且是无限重复的。 这意味着你的输入向量应该包含正弦信号的完整周期的整数个周期。您有64个采样的输入矢量和10的采样率。这会导致正弦波的6.4 周期,从而导致泄漏。如果您在执行FFT之后检查频谱 ,您会看到,没有两个干净的频率线,但两个地方附近有很多频率分量。
  2. 纠正你的输入矢量后,应该只有两个单独的频率值 不接近于零。 62个频率分量将由非常接近零的数字噪声组成。计算这些值的阶段会导致垃圾数据。
  3. 如果N是 输入采样的数量,则时域中pi/2的相移等效于时域的移位N/4。

我修改了你的代码。你会在下面找到它。 使用变量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; 

input signal three periods of a sine signal

输入信号与正弦信号的三个周期

spectrum of input signal. Frequency lines at -3 and +3

输入信号的频谱。在-3频率线和3

input signal (blue) and phase shifted signal (red)

输入信号(蓝色)和相移信号(红色)

+0

非常感谢!我非常感谢! – Kristian 2014-09-22 14:59:45