频域分析方法代码实现！前半段/后半段频谱截取怎样最好？

在信号处理领域，频域分析是至关重要的，但很多人在编程实现这一过程中遇到了困难。这恰恰是我们今天需要重点探讨的核心议题。

直接FFT的结果与调整

在执行FFT分析后，结果呈现出了明显的规律。例如，信号的起始部分对应的是频率区间[0,fs/2]，而结束部分则对应[-fs/2,0]。在实际操作中，必须将零频点放置在频谱的中央，这非得借助fftshift函数不可。以一个具体项目为例，在分析一段采集到的音频信号时，未经fftshift处理的FFT结果让人难以直接识别频率分布。但经过调整，频谱图便清晰地展现了正负频率的分布。这一步骤对于声音信号的频域分析等类似工作极为关键。

我们得知，处理FFT结果需要采用截取正频段的方法。这主要有两种方式：一是在fftshift后取后半部分，二是直接在fft结果中取前半部分。这两种方法得出的结果并无差异。就好比分析心电图信号，不管选用哪种截取方式，都是为了提取正频段成分，以便进行与疾病相关的频率特征分析。

t_s = 0.01; %采样周期
t_start = 0.5; %起始时间
t_end = 5;     %结束时间
t = t_start : t_s : t_end;
y = 1.5*sin(2*pi*5*t)+3*sin(2*pi*20*t)+randn(1,length(t));  %生成信号

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%频谱%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y_f = fft(y); %傅里叶变换
subplot(5,1,1);
plot(t,y);title('original signal');   %绘制原始信号图
Druation = t_end -t_start;  %计算采样时间
Sampling_points = Druation/t_s +1;  %采样点数，fft后的点数就是这个数
f_s = 1/t_s; %采样频率
f_x = 0:f_s/(Sampling_points -1):f_s;  %注意这里和横坐标频率对应上了，频率分辨率就是f_s/(Sampling_points -1)
t2 = f_x-f_s/2;
shift_f = abs(fftshift(y_f));
subplot(5,1,2);
plot(f_x,abs(y_f));title('fft transform');
subplot(5,1,3);
plot(f_x-f_s/2,shift_f);title('shift fft transform');  %将0频率分量移到坐标中心
subplot(5,1,4);
plot(t2(length(t2)/2:length(t2)),shift_f(length(shift_f)/2:length(shift_f)));title('shift fft transform');   %保留正频率部分
subplot(5,1,5);
plot(f_x(1:length(f_x)/2),abs(y_f(1:length(f_x)/2)));title('fft cut');  %直接截取fft结果的前半部分

奈奎斯特定理的作用

奈奎斯特定理在频域分析中扮演着核心角色。它指出，信号的采样速度必须高于其最高频率的两倍。比如，对于雷达信号，如果雷达发射的最高频率是10MHz，那么采样速度至少要是20MHz。如果低于这个标准，就会导致频谱重叠。实验表明，如果采样速度不够，采集到的频谱图中会出现错误频率，这会影响到所有频域处理结果的精确度。

换个角度审视，此定理明确了频域分析中基础采样准则的设定。比如，在处理地震波信号时，若不依照奈奎斯特定理，低频地震信号可能会与高频噪声混淆，进而造成地层结构关键信息的识别困难。

功率谱的求法

计算功率谱有俩种方法。一种是用傅立叶变换的平方除以区间长度来算，这方法挺简单，容易理解。在MATLAB编程里，我们根据这个公式就能得到结果。例如，针对描述粒子运动轨迹的离散信号，用这个方法算功率谱，就能看出粒子在不同频率上的能量分布情况。

还有一法是采用自相关函数执行傅里叶转换。从理论层面分析，这两种手段的成果本应相同。但实际操作中，例如在分析含噪电子设备输出信号时，我们发现自相关函数法在降噪效果上更为突出，绘制的曲线也更加平滑。这好比在喧闹场合努力辨识对方说话，运用自相关函数傅里叶转换这一技术，能更高效地去除噪音干扰，进而更清晰地呈现出信号的功率谱特征。

Fs = 1000;
nfft = 1000;  %fft采样点数

%产生序列
n = 0:1/Fs:1;
xn = cos(2*pi*100*n) + 3*cos(2*pi*200*n)+(randn(size(n)));
subplot(5,1,1);plot(xn);title('加噪信号');xlim([0 1000]);grid on
%FFT
Y = fft(xn,nfft);
Y = abs(Y);
subplot(5,1,2);plot((10*log10(Y(1:nfft/2))));title('FFT');xlim([0 500]);grid on
%FFT直接平方
Y2 = Y.^2/(nfft);
subplot(5,1,3);plot(10*log10(Y2(1:nfft/2)));title('直接法');xlim([0 500]);grid on
%周期图法
window = boxcar(length(xn));  %矩形窗
[psd1,f] = periodogram(xn,window,nfft,Fs);
psd1 = psd1 / max(psd1);
subplot(5,1,4);plot(f,10*log10(psd1));title('周期图法');ylim([-60 10]);grid on
%自相关结果
cxn = xcorr(xn,'unbiased');  %计算自相关函数
%自相关法
CXk = fft(cxn,nfft);
psd2 = abs(CXk);
index = 0:round(nfft/2-1);
k = index*Fs/nfft;
psd2 = psd2/max(psd2);
psd2 = 10*log10(psd2(index+1));
subplot(5,1,5);plot(k,psd2);title('间接法');grid on

直接法与周期图法对比

编写代码时，我们运用直接法，通过傅立叶变换的平方除以区间长度来求得功率谱。这一方法与MATLAB中的periodogram函数所得结果相同。在分析机械振动信号时，我们对比了这两种计算功率谱的方法。结果显示，在数据量不多且信号频率成分简单的情况下，两种方法得出的结果非常接近。这一发现为我们在编写频域分析代码时提供了更多选择，使我们能根据实际情况挑选出更易理解或计算速度更快的方案。

工业监测系统若要依赖实时功率谱来预测故障，这一点极为关键。以风力发电机的振动监控为例，我们需根据硬件条件和信号特点，巧妙选择计算功率谱的方法。最重要的是，要能快速且精确地识别出潜在的危险频率信号。

倒频谱的计算与含义

实倒频谱函数用于计算倒频谱，MATLAB手册中指出其计算步骤是先将信号转换为频谱，接着转为对数形式，最后执行傅立叶逆变换。但倒频谱的定义是将信号转成功率谱，随后进行对数处理，再进行傅立叶逆变换。两者间的差异可能源于功率谱是频谱值的平方，对数处理后，平方变成了系数2，但这对于后续计算影响微乎其微，所以近似结果并无差异。

在仿真实验里，我们需要留意倒频谱的作用，需要亲手制作一组调频信号。举例来说，对那些高频（如50Hz、100Hz、200Hz）和低频（如5Hz、10Hz、20Hz）的信号进行调整，接着分别画出低频、高频和调频信号的时域图和频谱图。这种调频信号的处理方法，在通信信号分析领域应用很广，有助于提高信号传输的效率和稳定性。

倒频谱在信号处理中的体现

sf = 1000;
nfft = 1000;
x = 0:1/sf:5;
y1=10*cos(2*pi*5*x)+7*cos(2*pi*10*x)+5*cos(2*pi*20*x)+0.5*randn(size(x));
y2=20*cos(2*pi*50*x)+15*cos(2*pi*100*x)+25*cos(2*pi*200*x)+0.5*randn(size(x));
for i = 1:length(x)
    y(i) = y1(i)*y2(i);
end
subplot(3,3,1)
plot(y1);xlim([0 5000]);title('y1');
subplot(3,3,2)
plot(y2);xlim([0 5000]);title('y2');
subplot(3,3,3)
plot(y);xlim([0 5000]);title('y=y1*y2');

t = 0:1/sf:(nfft-1)/sf;
nn = 1:nfft;
subplot(3,3,4)
ft = fft(y1,nfft);
Y = abs(ft);
plot(0:nfft/2-1,((Y(1:nfft/2))));
title('fft_y_1');
ylabel('幅值');xlim([0 300]);
grid on;
subplot(3,3,5)
ft = fft(y2,nfft);
Y = abs(ft);
plot(0:nfft/2-1,((Y(1:nfft/2))));
title('fft_y_2');
ylabel('幅值');xlim([0 300]);
grid on;
subplot(3,3,6)
ft = fft(y,nfft);
Y = abs(ft);
plot(0:nfft/2-1,((Y(1:nfft/2))));
title('fft_y');
ylabel('幅值');xlim([0 300]);
grid on;

subplot(3,3,7)
z = rceps(y);
plot(t(nn),abs(z(nn)));
title('z=rceps(y)');ylim([0 0.3]);
xlabel('时间(s)');
ylabel('幅值');
grid on;
subplot(3,3,8)
yy = real(ifft(log(abs(fft(y))))); %信号→傅里叶→对数→傅里叶逆变换
plot(t(nn),abs(yy(nn)));
title('real(ifft(log(abs(fft(y)))))');ylim([0 0.3]);
xlabel('时间(s)');
ylabel('幅值');
grid on;

以特定调制信号为例，例如某些信号中带有20Hz、10Hz和5Hz等低频部分。在常规的FFT_y分析中，这些低频部分以边缘频带的形式出现，难以准确判断其频率。但在倒频谱分析中，这些部分却很容易被识别。这好比在黑夜中，一双隐藏的眼睛在FFT_y分析中难以发现，而在倒频谱分析中却能立刻被“照亮”，显示出其特征。在无线电通信的信号分析中，倒频谱分析这种处理类似问题的能力对工程师来说非常宝贵，它能够帮助他们快速发现并解读微弱的低频信号。

关于项目操作中选用频域分析方法的问题，我想听听大家的意见：面对不同项目的具体特点，你们是如何做出选择的？欢迎在评论区留下你们的想法。若这篇文章对您有所启发，不妨点赞或分享给更多人。