频域特征提取的Python实现（频谱、功率谱、倒频谱）MATLAB程序代码：
%==========================================================================
%Desc:      以⾼斯信号为例，求解其频谱、双边功率谱、单边功率谱、双边功率谱密度、
%          单边功率谱密度，这⾥⾼斯信号的半波全宽FWHM=50ps，中⼼点位于2.5ns处。
%=========================================================================
clc;
clear;
FWHM=50e-12;            %⾼斯信号FWHM宽度，为50ps
time_window=100*FWHM;  %⾼斯信号的采样窗⼝宽度，该值决定了傅⾥叶变换后的频率分辨率
Ns=2048;                %采样点
dt=time_window/(Ns-1);  %采样时间间隔
t=0:dt:time_window;    %采样时间
gauss_time=exp(-0.5*(2*sqrt(2*log(2))*(t-2.5e-9)/FWHM).^2); %⾼斯脉冲，中⼼位于2.5ns处。
plot(t*1e+9,gauss_time,'linewidth',2.5);
xlabel('Time/ns');
ylabel('Amplitude/V');
title('Gauss pulse');
%===========以下计算双边谱、双边功率谱、双边功率谱密度=================
gauss_spec=fftshift(fft(ifftshift(gauss_time)));    %傅⾥叶变换，并且进⾏fftshift移位操作。
gauss_spec=gauss_spec/Ns;  %求实际的幅度值；
df=1/time_window;              %频率分辨率
k=floor(-(Ns-1)/2:(Ns-1)/2);
% k=0:Ns-1;
double_f=k*df;  %双边频谱对应的频点
figure; %幅度谱
plot(double_f*1e-9,abs(gauss_spec),'linewidth',2.5);
xlabel('Frequency/GHz');
ylabel('Amplitude/V');
title('double Amplitude spectrum');
figure; %相位谱
plot(double_f*1e-9,angle(gauss_spec),'linewidth',2.5);
xlabel('Frequency/GHz');
ylabel('Phase/rad');
title('double Phase spectrum');
figure; %功率谱
double_power_spec_W=abs(gauss_spec).^2;                %双边功率谱，单位W；
double_power_spec_mW=double_power_spec_W*1e+3;          %双边功率谱，单位mW；
double_power_spec_dBm=10*log10(double_power_spec_mW);  %双边功率谱，单位dBm；
plot(double_f*1e-9,double_power_spec_dBm,'linewidth',2.5);
xlabel('Frequency/GHz');
ylabel('Power/dBm');
title('double Power spectrum');
figure; %功率谱密度
double_power_specD_W=abs(gauss_spec).^2/(df);      %双边功率谱密度，单位W/Hz
double_power_specD_mW=double_power_specD_W*1e+3;    %双边功率谱密度，单位mW/Hz
double_power_specD_dBm=10*log10(double_power_specD_mW);%双边功率谱密度，单位dBm/Hz
plot(double_f*1e-9,double_power_specD_dBm,'linewidth',2.5);
xlabel('Frequency/GHz');
ylabel('Power/(dBm/Hz)');
title('double power spectrum Density');
%==========以下计算单边谱、单边功率谱及单边功率谱密度=========
%==========以下计算单边谱、单边功率谱及单边功率谱密度=========
gauss_spec=fft(ifftshift(gauss_time));  %计算单边谱⽆需fftshift
gauss_spec=gauss_spec/Ns;      %计算真实的幅度值
matlab傅里叶变换的幅度谱和相位谱single_gauss_spec=gauss_spec(1:floor(Ns/2));
single_f=(0:floor(Ns/2)-1)*df;
figure; %幅度谱
plot(single_f*1e-9,abs(single_gauss_spec),'linewidth',2.5);
xlabel('Frequency/GHz');
ylabel('Amplitude/V');
title('single Amplitude spectrum');
figure; %相位谱
plot(single_f*1e-9,angle(single_gauss_spec),'linewidth',2.5);
xlabel('Frequency/GHz');
ylabel('Phase/rad');
title('single Phase spectrum');
figure;%功率谱
double_power_spec_W=abs(gauss_spec).^2;
single_power_spec_W=2*double_power_spec_W(1:floor(Ns/2));  %单边功率谱，单位W
single_power_spec_mW=single_power_spec_W*1e+3;              %单边功率谱，单位mW；
single_power_spec_dBm=10*log10(single_power_spec_mW);      %双边功率谱，单位dBm；
plot(single_f*1e-9,single_power_spec_dBm,'linewidth',2.5);
xlabel('Frequency/GHz');
ylabel('Power/dBm');
title('single Power spectrum');
figure;%功率谱密度
double_power_specD_W=abs(gauss_spec).^2/(df);
single_power_specD_W=2*double_power_specD_W(1:floor(Ns/2));        %单边功率谱密度，单位W/Hz single_power_specD_mW=single_power_specD_W*1e+3;                    %单边功率谱密度，单位mW/Hz single_power_specD_dBm=10*log10(single_power_specD_mW);            %单边功率谱密度，单位dBm/Hz plot(single_f*1e-9,single_power_specD_mW,'linewidth',2.5);
xlabel('Frequency/GHz');
ylabel('Power/(dBm/Hz)');
title('single power spectrum density');
python 代码
频谱
from scipy.fftpack import fft, fftshift, ifft
from scipy.fftpack import fftfreq
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
t_s = 0.01
t_start = 0.5
t_end = 5
t = np.arange(t_start, t_end, t_s)
f0 = 5
f1 = 20
# generate the orignal signal
y = 1.5*np.sin(2*np.pi*f0*t) + 3*np.sin(2*np.pi*20*t) + np.random.randn(t.size)
# fft
Y = fft(y)
# fftshift
shift_Y = fftshift(Y)
# the positive part of fft, get from fft
pos_Y_from_fft = Y[:Y.size//2]
# the positive part of fft, get from shift fft
pos_Y_from_shift = shift_Y[shift_Y.size//2:]
# plot the figures
plt.figure(figsize=(10, 12))
plt.subplot(511)
plt.plot(y)
plt.subplot(512)
plt.plot(np.abs(Y))
plt.subplot(513)
plt.plot(np.abs(shift_Y))
plt.subplot(514)
plt.plot(np.abs(pos_Y_from_fft))
plt.subplot(515)
plt.plot(np.abs(pos_Y_from_shift))
plt.show()
功率谱
# generate original signal
t = np.arange(0, 1, 1/fs)
f0 = 100
f1 = 200
x = np.cos(2*np.pi*f0*t) + s(2*np.pi*f1*t) + np.random.randn(t.size)
# FFT
Y = fft(x, num_fft)
Y = np.abs(Y)
# power spectrum
ps = Y**2 / num_fft
# power spectrum using correlate
cor_x = np.correlate(x, x, 'same')
cor_X = fft(cor_x, num_fft)
ps_cor = np.abs(cor_X)
ps_cor = ps_cor / np.max(ps_cor)
plt.figure(figsize=(15, 12))
plt.subplot(511)
plt.plot(x)
plt.subplot(512)
plt.plot(20*np.log10(Y[:num_fft//2]))
plt.subplot(513)
plt.plot(20*np.log10(ps[:num_fft//2]))
plt.subplot(514)
plt.plot(20*np.log10(ps_cor[:num_fft//2]))
plt.show()
倒频谱
t = np.arange(0, 5, 1/fs)
y1 = s(2*np.pi*5*t) + s(2*np.pi*10*t) + s(2*np.pi*20*t) + np.random.randn(t.size)
y2 = s(2*np.pi*50*t) + s(2*np.pi*100*t) + s(2*np.pi*200*t) + np.random.randn(t.size) y = y1*y2
Y1 = fft(y1, num_fft)
Y1 = np.abs(Y1)
Y2 = fft(y2, num_fft)
Y2 = np.abs(Y2)
Y = fft(y, num_fft)
Y = np.abs(Y)
spectrum = np.fft.fft(y, n=num_fft)
ceps = np.fft.ifft(np.log(np.abs(spectrum))).real
plt.figure(figsize=(20, 12))
plt.subplot(331)
plt.plot(y1)
plt.subplot(332)
plt.plot(y2)
plt.subplot(333)
plt.plot(y)
plt.subplot(334)
plt.plot(Y1[:num_fft//2])
plt.subplot(335)
plt.plot(Y2[:num_fft//2])
plt.subplot(336)
plt.plot(Y[:num_fft//2])
plt.show()

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。