频谱分析

计算频谱和功率谱密度

function main
close all;
clear all;
clc;

data = load('D:Desktop\datastate1data212.txt');
t = load('D:Desktop\data	ime.txt');

Ts = 3.1e-3;
Fs = 1/Ts;
% plot(t, data(:, 5)*9.8)

ch1 = data(:, 1); % 单位N
ch2 = data(:, 2); % 单位g
ch3 = data(:, 3); % 单位g
ch4 = data(:, 4); % 单位g
ch5 = data(:, 5); % 单位g

% 加窗
% wvtool(hann(128), hamming(128), gausswin(128)) % 窗可视化工具
% wintool % 窗设计工具

win = hanning(length(t));
ch1win = ch1.*win;
ch2win = ch2.*win*9.8;
ch3win = ch3.*win*9.8;
ch4win = ch4.*win*9.8;
ch5win = ch5.*win*9.8;

% 不加窗
[Yf, f] = Spectrum_Calc(ch5, Fs);
figure
subplot(2, 1, 1)
plot(t, ch5)
grid on
xlabel('t')
ylabel('channel-1')
title('原始信号')
subplot(2, 1, 2)
plot(f, Yf)
grid on
xlabel('f')
ylabel('|Yf|')
% xlim([0 100])
% ylim([0 1])
title('原始信号频谱')

% 加窗
[Yf_win, f_win] = Spectrum_Calc(ch5win, Fs);
figure
subplot(2, 1, 1)
plot(t, ch5win)
grid on
xlabel('t')
ylabel('y')
title('加窗信号')
subplot(2, 1, 2)
plot(f_win, 2*Yf_win) % 2表示能量系数
grid on
xlabel('f')
ylabel('|Yf|')
% xlim([0 100])
% ylim([0 1])
title('加窗信号频谱')

%% welch法计算功率谱密度（PSD）
N = 1024;
figure
pwelch(ch2, hanning(N),N/2, N, Fs) % 分段加窗，重叠1/2
end
%% 求取频谱
function [Yf, f] = Spectrum_Calc(yt, Fs)
    L = length(yt);
    NFFT = 2^nextpow2(L);
    Yf = fft(yt,NFFT)/L;
    Yf = 2*abs(Yf(1:NFFT/2+1));
    f = Fs/2*linspace(0,1,NFFT/2+1);
end