matlab设计复合信号不同频率的初相位

生成原始信号

1. 为了方便起见，生成一个包含三个频率的复信号，分别是5Hz、10Hz、15Hz，初相位分别为采样率取100，采样点取512个，代码如下。

   clc;clear;
   Fs =100;%采样率
   N=512;%序列长度
   T = 1/Fs;%采样间隔
   t = 0:T:(N-1)*T;%时间序列
   s1 = cos(2*pi*5.*t+pi/4)+cos(2*pi*10.*t+pi/8)+cos(2*pi*15.*t+pi/4);%信号实部
   s2 = sin(2*pi*5.*t+pi/4)+sin(2*pi*10.*t+pi/8)+sin(2*pi*15.*t+pi/4);%信号虚部
   ss = complex(s1,s2);%合成复信号
   

使用FFT的方法

1. 对信号做fft，并生成相位谱和幅度谱：

   y = fft(ss);
   f = Fs.*((0:n-1)-(n/2))./(n);
   ys = fftshift(y);
   plot(f,abs(ys));%画出幅度谱
   

相位谱如下所示：

image

但是发现5，10，15赫兹的位置幅值不同。再来算一下15Hz处的相位：

pos1 = find(abs(y)==max(abs(y)));%这里利用了15Hz时幅值最大
tempp = angle(y);
tempp3 = tempp(pos1);

输出结果为0.1497弧度，显然不等于。

可能的原因：应该时由于信号长度取为512，它不是信号的周期。从而导致信号在周期性延拓后变得与原信号不一样。修改N的值为600(信号的周期)，重复刚才的分析：

clc;clear;
Fs =100;
N=600;%注意这里
T = 1/Fs;
t = 0:T:(N-1)*T;
s1 = cos(2*pi*5.*t+pi/4)+cos(2*pi*10.*t+pi/8)+cos(2*pi*15.*t+pi/4);
s2 = sin(2*pi*5.*t+pi/4)+sin(2*pi*10.*t+pi/8)+sin(2*pi*15.*t+pi/4);
ss = complex(s1,s2);
y = fft(ss);
f = Fs.*((0:n-1)-(n/2))./(n);
ys = fftshift(y);
plot(f,abs(ys));%画出幅度谱
pos1 = find(abs(y)==max(abs(y)));%此时对应的是10Hz
tempp = angle(y);
tempp3 = tempp(pos1);

幅度谱为：

image

​ 可以看到此时幅度谱正常，10Hz处的相位计算结果为0.3927，为预期值。

2. 总结：

   • 以信号的周期长度截取信号

   • 对信号做fft

   • 计算待测频率对应的横坐标，设是待求的频点，则对应fft序列中的第个数，满足：
     
     其中是采样率，是fft序列长度。

使用fir设计梳状滤波器，滤除其他频率

1. 设计梳妆滤波器

   Fo = [10 15];%待滤除的频率
   Wo = Fo./(Fs/2);%将频率归一化
   Wo = [0 Wo 0.3 0.5 1]%添加一些其他的频点要求，使得序列长度为偶数
   %生成一个与Wo等长的矩阵，滤波器的设计目的是使得滤波后的信号，在Wo的频点处的增益情况趋近于aa中对应的值
   %在这里让Fo对应的值为0，其他为1。
   aa = ones(1,length(Wo));
   aa(1,2)=0;
   aa(1,3)=0;
   %生成一个与aa等长的矩阵，对aa处对应的值进行约束。这里如果aa为1，则取'n'，为0则取's'。详情参见'help firgr'
   for i = 1:length(aa)
       if(aa(1,i)==1)
           s = [s 'n'];
       else
           s = [s 's']
       end
   end
   b = firgr(30,Wo,aa,s);%生成滤波器，第一个参数是阶数的要求，需要是个偶数
   freqz(b,1);%查看滤波器的相位谱和幅度谱
   sf = filter(b,1,ss);%对信号进行滤波
   

   频谱图如下：

   image

   滤波后的信号幅度谱如下：

   image

   可以看到5Hz和10Hz的频率基本没有了。

   但是我们如果直接使用angle(sg)来查看信号的相位信息会不准，原因是滤波器在滤波的时候添加了一个相位偏移，这点在滤波器的相位谱中可以体现出来。因此计算信号的初相位需要将这个相位偏移给减掉。

   使用命令[h1,h2]=freqz(b,1);来获得滤波器的频域响应，其中h1是频域的函数，angle(h1)可以获得对应频点上的频偏。

2. 计算频率在信号fft之后对应第几个点.

   假设做点fft，则对横坐标对应的频率为，其中是采样率。因此如果待求频率为，则它对应的fft的横坐标为（注意matlab下标从1开始）：
   
   其中route表示四舍五入。