信号频谱分析 功率谱分析 倒谱分析 小波分析

PDF获取路径

1.频谱分析

在本科信号系统课程中学习过傅里叶变换，可将信号时域波形转化为频域。为什么要进行域转换呢？因为大部分信号在传输过程中可能会受到外界因素的干扰（可以理解为"噪声"），这种干扰在时域上表现得不太明显，因此可以通过傅立叶变换将原来难以处理的时域信号转换成了易于分析的频域信号（信号的频谱）。

傅立叶原理表明：任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加。而根据该原理创立的傅立叶变换算法利用直接测量到的原始信号，以累加方式来计算该信号中不同正弦波信号的频率、振幅和相位。和傅立叶变换算法对应的是反傅立叶变换算法。该反变换从本质上说也是一种累加处理，这样就可以将单独改变的正弦波信号转换成一个信号。

clear; clc;t_s = 0.01; % 采样周期, 采样频率fs为100Hzt_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); % 傅里叶变换figure;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);len_t2 = floor(length(t2)/2);len_f = floor(length(shift_f)/2);plot(t2(len_t2:len_t2*2),shift_f(len_f:len_f*2));title('shift fft transform'); % 保留正频率部分subplot(5,1,5);len_fx = floor(length(f_x)/2);plot(f_x(1:len_fx),abs(y_f(1:len_fx)));title('fft cut'); % 直接截取fft结果的前半部分

第一幅图为生成信号的时域波形图，频率成分为5、20Hz，含一定程度噪声干扰。

第二幅图为该信号经傅里叶变换得到的频谱图，中心轴对称；

第三幅图为图二频谱搬移情况，将0频率分量搬到中心；

第四幅图为图三正频率部分；

第五幅图为图二结果的前半部分。

2. 功率谱分析

另外，可以从含噪信号功率谱中观察各个子频率强度；求取功率谱有fft法、周期图法和自相关法三种方式。

% 1.(傅立叶变换的平方)/(区间长度)；->直接法% 2.自相关函数的傅里叶变换；->相关函数法clear; clc;Fs = 1000; % 采样频率nfft = 1000; % 采样点数% 产生序列n = 0:1/Fs:1;xn = cos(2*pi*100*n) + 3*cos(2*pi*200*n) + (randn(size(n)));figure;subplot(5,1,1); plot(xn);title('加噪信号'); xlim([0,1000]); grid on;% FFTY = abs(fft(xn,nfft));subplot(5,1,2); plot((10*log10(Y(1:nfft/2))));title('FFT'); xlim([0,500]); grid on;% 法1：FFT直接平方Y2 = Y.^2/(nfft);subplot(5,1,3); plot(10*log10(Y2(1:nfft/2)));title('直接法'); xlim([0,500]); grid on;% 法2：周期图法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;% 法3：自相关结果cxn = xcorr(xn,'unbiased'); % 计算自相关函数% 法4：自相关法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;% 取对数的目的是使那些振幅较低的成分相对高振幅成分得以拉高，% 以便观察掩盖在低幅噪声中的信号特征。

从图中看出通过自相关法得到的功率谱成分（100Hz，200Hz）更加明显，其它频带干扰较少。

3. 倒谱分析

科学上网

blog

clear; clc;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); % 被调制endfigure;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)Y = abs(fft(y1,nfft)); % 1024点fftplot(0:nfft/2-1,((Y(1:nfft/2))));title('fft_y_1');ylabel('幅值');xlim([0 300]);grid on;subplot(3,3,5)Y = abs(fft(y2,nfft));plot(0:nfft/2-1,((Y(1:nfft/2))));title('fft_y_2');ylabel('幅值');xlim([0 300]);grid on;subplot(3,3,6)Y = abs(fft(y,nfft));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;

可以发现第7幅图跟第8幅图一样，说明了倒谱等价于信号→傅里叶→对数→傅里叶逆变换。

进阶

举一个例子：

通过倒谱达到语音分离效果；

注意：这里的时间并不是严格意义上的时间，而是伪频率。

load mtlb% soundsc(mtlb,Fs);% To hear, type soundsc(mtlb,Fs)timelag = 0.23; % 时延delta = round(Fs*timelag);alpha = 0.5; % 衰减系数orig = [mtlb;zeros(delta,1)];echo = [zeros(delta,1);mtlb]*alpha;mtEcho = orig + echo;pause(0.5);soundsc(mtEcho,Fs);t = (0:length(mtEcho)-1)/Fs;figure;subplot(2,1,1)plot(t,[orig echo])legend('Original','Echo')subplot(2,1,2)plot(t,mtEcho)legend('Total')xlabel('Time (s)')c = rceps(mtEcho);[px,locs] = findpeaks(c,'Threshold',0.2,'MinPeakDistance',0.2);clfplot(t,c,t(locs),px,'o')xlabel('Time (s)')dl = locs(2)-1;mtNew = filter(1,[1 zeros(1,dl-1) alpha],mtEcho);subplot(2,1,1)plot(t,orig)legend('Original')subplot(2,1,2)plot(t,mtNew)legend('Filtered')xlabel('Time (s)')% To hear, type soundsc(mtNew,Fs)soundsc(mtNew,Fs)

于是，两个语音段得到了分离；

4. 小波分析

clear; clc;t_s = 0.005; % 采样周期t_start = 0.001; % 起始时间t_end = 10; % 结束时间t = t_start : t_s : t_end;y = 10*sin(2*pi*2*t)+3*sin(2*pi*10*t)+1*sin(2*pi*20*t)+3*randn(1,length(t)); %生成信号len = length(y);% 生成突变信号y2 = 50*sin(2*pi*50*t); % 高振幅高频率 突变信号for i = 1: lenif i>=601&&i<=604y(i) = y(i)+y2(i);elsey(i) = y(i);endendfigure;plot(y) % 绘制原始信号[c,l] = wavedec(y,5,'db5');% 重构1~5层细节函数d5 = wrcoef('d',c,l,'db5',5);d4 = wrcoef('d',c,l,'db5',4);d3 = wrcoef('d',c,l,'db5',3);d2 = wrcoef('d',c,l,'db5',2);d1 = wrcoef('d',c,l,'db5',1);% 重构1~5层近似函数a5 = wrcoef('a',c,l,'db5',5);a4 = wrcoef('a',c,l,'db5',4);a3 = wrcoef('a',c,l,'db5',3);a2 = wrcoef('a',c,l,'db5',2);a1 = wrcoef('a',c,l,'db5',1);figuresubplot(5,2,1);plot(a1)subplot(5,2,2);plot(d1)subplot(5,2,3);plot(a2)subplot(5,2,4);plot(d2)subplot(5,2,5);plot(a3)subplot(5,2,6);plot(d3)subplot(5,2,7);plot(a4)subplot(5,2,8);plot(d4)subplot(5,2,9);plot(a5)subplot(5,2,10);plot(d5)figure;shift_fa5 = abs(fft(a5,2000)); % 傅里叶变换shift_fd5 = abs(fft(d5,2000)); % 傅里叶变换shift_fd4 = abs(fft(d4,2000)); % 傅里叶变换shift_fd3 = abs(fft(d3,2000)); % 傅里叶变换shift_fd2 = abs(fft(d2,2000)); % 傅里叶变换shift_fd1 = abs(fft(d1,2000)); % 傅里叶变换plot(shift_fa5(2:1001),'r');title('shift fft transform');hold on;plot(shift_fd5(2:1001),'g');title('shift fft transform');hold on;plot(shift_fd4(2:1001),'b');title('shift fft transform');hold on;plot(shift_fd3(2:1001),'y');title('shift fft transform');hold on;plot(shift_fd2(2:1001),'r');title('shift fft transform');hold on;plot(shift_fd1(2:1001),'g');title('shift fft transform');hold on;

小波分析是一种基于频率逐步二分的方法，信号x=a1+d1；a1、d1分量频率成分各为一半；a1=a2+d2，a2、d2分量频率成分各为一半，以此类推。

详细了解请科学上网

小波分解得到各成分的频率分布为：

5. 时域特征值

% 均方根、峰值因子、脉冲因子、裕度因子、峭度因子、波形因子和偏度clear; clc;% x=0:0.1:2*pi;% y=sin(x);%信号y = normrnd(0,2,2000,1); % 生成平均值为0，标准差为2的1000*1的矩阵figure;plot(y);ma = max(y); % 最大值mi = min(y); % 最小值me = mean(y); % 平均值pk = ma-mi; % 峰-峰值av = mean(abs(y));% 绝对值的平均值(整流平均值)va = var(y);% 方差st = std(y);% 标准差ku = kurtosis(y);% 峭度sk = skewness(y); % 偏度rm = rms(y);% 均方根S = rm/av; % 波形因子C = pk/rm; % 峰值因子I = pk/av; % 脉冲因子xr = mean(sqrt(abs(y)))^2;L = pk/xr; % 裕度因子