生成原始信号

为了方便起见，生成一个包含三个频率的复信号，分别是5Hz、10Hz、15Hz，初相位分别为 $\frac{\pi}{4},\frac{\pi}{8},\frac{\pi}{4}$ 采样率取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的方法

对信号做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弧度，显然不等于 $\frac{\pi}{4}$ 。

可能的原因：应该时由于信号长度取为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，为预期值 $\frac{pi}{8}$ 。

总结：
- 以信号的周期长度截取信号
- 对信号做fft
- 计算待测频率对应的横坐标，设 $f$ 是待求的频点，则对应fft序列中的第 $i$ 个数，满足：
  $\frac{f}{Fs}=\frac{i-1}{N}$
  其中 $Fs$ 是采样率， $N$ 是fft序列长度。

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

设计梳妆滤波器

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)可以获得对应频点上的频偏。

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

假设做 $n$ 点fft，则对横坐标对应的频率为 $(0:(n-1))./n.*Fs$ ，其中 $Fs$ 是采样率。因此如果待求频率为 $f$ ，则它对应的fft的横坐标为（注意matlab下标从1开始）：
$route(\frac{f*n}{Fs})+1$
其中route表示四舍五入。