生成原始信號(hào)

1. 為了方便起見，生成一個(gè)包含三個(gè)頻率的復(fù)信號(hào)，分別是5Hz、10Hz、15Hz，初相位分別為采樣率取100，采樣點(diǎn)取512個(gè)，代碼如下。

   clc;clear;
   Fs =100;%采樣率
   N=512;%序列長(zhǎng)度
   T = 1/Fs;%采樣間隔
   t = 0:T:(N-1)*T;%時(shí)間序列
   s1 = cos(2*pi*5.*t+pi/4)+cos(2*pi*10.*t+pi/8)+cos(2*pi*15.*t+pi/4);%信號(hào)實(shí)部
   s2 = sin(2*pi*5.*t+pi/4)+sin(2*pi*10.*t+pi/8)+sin(2*pi*15.*t+pi/4);%信號(hào)虛部
   ss = complex(s1,s2);%合成復(fù)信號(hào)
   

使用FFT的方法

1. 對(duì)信號(hào)做fft，并生成相位譜和幅度譜：

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

相位譜如下所示：

image

但是發(fā)現(xiàn)5，10，15赫茲的位置幅值不同。再來(lái)算一下15Hz處的相位：

pos1 = find(abs(y)==max(abs(y)));%這里利用了15Hz時(shí)幅值最大
tempp = angle(y);
tempp3 = tempp(pos1);

輸出結(jié)果為0.1497弧度，顯然不等于。

可能的原因：應(yīng)該時(shí)由于信號(hào)長(zhǎng)度取為512，它不是信號(hào)的周期。從而導(dǎo)致信號(hào)在周期性延拓后變得與原信號(hào)不一樣。修改N的值為600(信號(hào)的周期)，重復(fù)剛才的分析：

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)));%此時(shí)對(duì)應(yīng)的是10Hz
tempp = angle(y);
tempp3 = tempp(pos1);

幅度譜為：

image

? 可以看到此時(shí)幅度譜正常，10Hz處的相位計(jì)算結(jié)果為0.3927，為預(yù)期值。

2. 總結(jié)：

   • 以信號(hào)的周期長(zhǎng)度截取信號(hào)

   • 對(duì)信號(hào)做fft

   • 計(jì)算待測(cè)頻率對(duì)應(yīng)的橫坐標(biāo)，設(shè)是待求的頻點(diǎn)，則對(duì)應(yīng)fft序列中的第個(gè)數(shù)，滿足：
     
     其中是采樣率，是fft序列長(zhǎng)度。

使用fir設(shè)計(jì)梳狀濾波器，濾除其他頻率

1. 設(shè)計(jì)梳妝濾波器

   Fo = [10 15];%待濾除的頻率
   Wo = Fo./(Fs/2);%將頻率歸一化
   Wo = [0 Wo 0.3 0.5 1]%添加一些其他的頻點(diǎn)要求，使得序列長(zhǎng)度為偶數(shù)
   %生成一個(gè)與Wo等長(zhǎng)的矩陣，濾波器的設(shè)計(jì)目的是使得濾波后的信號(hào)，在Wo的頻點(diǎn)處的增益情況趨近于aa中對(duì)應(yīng)的值
   %在這里讓Fo對(duì)應(yīng)的值為0，其他為1。
   aa = ones(1,length(Wo));
   aa(1,2)=0;
   aa(1,3)=0;
   %生成一個(gè)與aa等長(zhǎng)的矩陣，對(duì)aa處對(duì)應(yīng)的值進(jìn)行約束。這里如果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);%生成濾波器，第一個(gè)參數(shù)是階數(shù)的要求，需要是個(gè)偶數(shù)
   freqz(b,1);%查看濾波器的相位譜和幅度譜
   sf = filter(b,1,ss);%對(duì)信號(hào)進(jìn)行濾波
   

   頻譜圖如下：

   image

   濾波后的信號(hào)幅度譜如下：

   image

   可以看到5Hz和10Hz的頻率基本沒(méi)有了。

   但是我們?nèi)绻苯邮褂?code>angle(sg)來(lái)查看信號(hào)的相位信息會(huì)不準(zhǔn)，原因是濾波器在濾波的時(shí)候添加了一個(gè)相位偏移，這點(diǎn)在濾波器的相位譜中可以體現(xiàn)出來(lái)。因此計(jì)算信號(hào)的初相位需要將這個(gè)相位偏移給減掉。

   使用命令[h1,h2]=freqz(b,1);來(lái)獲得濾波器的頻域響應(yīng)，其中h1是頻域的函數(shù)，angle(h1)可以獲得對(duì)應(yīng)頻點(diǎn)上的頻偏。

2. 計(jì)算頻率在信號(hào)fft之后對(duì)應(yīng)第幾個(gè)點(diǎn).

   假設(shè)做點(diǎn)fft，則對(duì)橫坐標(biāo)對(duì)應(yīng)的頻率為，其中是采樣率。因此如果待求頻率為，則它對(duì)應(yīng)的fft的橫坐標(biāo)為（注意matlab下標(biāo)從1開始）：
   
   其中route表示四舍五入。