频谱选带与细化分析两个小例子

频谱选带与细化分析两个小例子

解：

利用matlab编程，得出图9-1所示结果

（a）信号采样图

（b）参与FFT为64点频谱图

（c）参与FFT为128点频谱图

（d）参与FFT为256点频谱图

（e）参与FFT为64点补零至256点频谱图

（f）参与FFT为128点补零至256点频谱图

图9-1 第9章第1题各结果图示

作出一个周期内各幅值谱图如图9-1所示。

比较图（b）参与FFT为64个点，（c）参与FFT为128个点，（d）参与FFT为256个点。图（b）只有一个峰值，峰较宽；图（c）也只有一个峰值，峰较窄；图（d）有三个峰值，峰更窄。显然只有图（d）能分辨出信号中有三个频率（49Hz，50Hz，51Hz）。容易得出，随着参与点数增加，分辨率越来越高，这个是显然的。

比较图（b）参与FFT为64点频谱图，图（e）参与FFT为64点补零至256点频谱图，可以发现，补零之后的频谱分辨率提高了，但造成谱形失真。

比较图（c）参与FFT为128点频谱图和图（f）参与FFT为128点补零至256点频谱图。可以发现，补零之后的频谱分辨率提高了，能大致分辨出三个峰值，但造成谱形失真。

比较图（e）参与FFT为64点补零至256点频谱图和图（f）参与FFT为128点补零至256点频谱图。图（e）失真更严重，因此补零越多，失真越厉害。

该题的Matlab代码如下：

clear;clc;

%信号采样，采样频率128Hz，采样点数256

Fs=128;

N=256;

n=0:N-1;

t=n/Fs;

x=100*cos(2*pi*50*t)+150*cos(2*pi*51*t)+50*cos(2*pi*49*t);

subplot(3,2,1);

plot(t,x) %画出采样信号

xlabel('t/s'); ylabel('Amplitude'); title('信号采样图');

%参与FFT点数为64

N1=64;

y1=fft(x(1:N1));%对信号进行快速Fourier变换

amp1=abs(y1)*2/N1;%求Fourier变换振幅

f=n(1:N1)*Fs/N1;

subplot(3,2,3);

plot(f(1:N1/2),amp1(1:N1/2))

xlabel('f/Hz'); ylabel('Amplitude1'); title('参与FFT点64频谱图'); %参与FFT点数为128

N2=128;

y2=fft(x(1:N2));

amp2=abs(y2)*2/N2;

f=n(1:N2)*Fs/N2;

subplot(3,2,5);

plot(f(1:N2/2),amp2(1:N2/2))

xlabel('f/Hz'); ylabel('Amplitude2'); title('参与FFT点128频谱图'); %参与FFT点数为256

y3=fft(x(1:N));%对信号进行快速Fourier变换

amp3=abs(y3)*2/N;%求Fourier变换振幅

f=n(1:256)*Fs/N;

subplot(3,2,2);

plot(f(1:N/2),amp3(1:N/2))

xlabel('f/Hz'); ylabel('Amplitude3'); title('参与FFT点256频谱图'); %对64个点样本补零增为256点后，进行FFT变换

x1=[x(1:64),zeros(1,192)];

y4=fft(x1);

amp4=abs(y4)*2/N1;

f=n(1:256)*Fs/N;

subplot(3,2,4);

plot(f(1:N/2),amp4(1:N/2))

xlabel('f/Hz'); ylabel('Amplitude4'); title('参与FFT点64补零至256频谱图') %对128个点样本补零增为256点后，进行FFT变换

x1=[x(1:128),zeros(1,128)];

y5=fft(x1);

amp5=abs(y5)*2/N2;

f=n*Fs/N;

subplot(3,2,6);

plot(f(1:N/2),amp5(1:N/2))

xlabel('f/Hz'); ylabel('Amplitude5'); title('参与FFT点128补零至256频谱图');

解：

利用matlab编程，得出图9-2所示结果

（a）信号采集图

（b）全景频谱图

（c）细化后频谱图

图9-2 第9章2题结果图示

表9-1 计算结果

参数理论值FFT谱细化谱频率（Hz）84.6 84.0000 84.6000 幅值100 85.7550 99.8996 相位（度）54 107.1010 65.4522

观察图9-2易知细化后的谱图更加清楚。通过表9-1中的计算结果，易知，细化谱的频率（Hz）、幅值、相位（度）与理论值更加接近，细化的目的达到。

该题Matlab代码如下：

clear;clc;

%信号采样，采样频率Fs=256Hz，采样点数N=128

Fs=256;

N=128;

n=0:N-1;

t=n/Fs;

x=100*cos(2*pi*84.6*t+0.3*pi);

%画信号采样图

subplot(3,1,1);

plot(t,x)

xlabel('t/s'); ylabel('Amplitude1'); title('信号采样图');

%全景FFT变换

y1=fft(x);%对信号进行快速Fourier变换

amp1=abs(y1).*2./N;%求Fourier变换振幅（需要将结果的绝对值乘以2除以N）Theta1=angle(y1)*180/pi;%计算相位

% Theta1=atan(imag(y1)./real(y1));

kmax1=find(amp1==max(amp1)); %找出幅值最大的序列下标值

f=n*Fs/N;%归一化频率轴

%作全景幅值谱图

subplot(3,1,2);

plot(f(1:N/2),amp1(1:N/2))

xlabel('f/Hz'); ylabel('Amplitude2');title('全景频谱图');

fprintf(' 全景幅值谱图最大幅值为:%3.4f\n ',amp1(kmax1(1)))

fprintf('全景幅值谱图最大幅值处频率为:%3.4fHz\n ',f(kmax1(1)))

fprintf('全景幅值谱图最大幅值处相位为:%3.4f度\n ',Theta1(kmax1(1)));

%对82~86Hz之间进行频谱细化，取谱线条数为41条

fx=82:0.1:86;%取41条频率

for i=1:length(fx)

af(i)=2/N*x*cos(2*pi*(1:N)*fx(i)/Fs)';%计算a(f)的值，即实部

bf(i)=2/N*x*sin(2*pi*(1:N)*fx(i)/Fs)';%计算b(f)的值，即虚部

amp2(i)=sqrt(af(i)^2+bf(i)^2);%计算幅值

Theta2(i)=angle(af(i)+bf(i)*1i)*180/pi;%计算相位

% Theta2(i)=atan((-bf(i))/af(i));

end