clear, clc, close all
%% example1
% 音频信号处理实例
% [x, fs] = audioread('track.wav'); % 微软开机声音
% x = x(:, 1);
% xlen = length(x);
%
% t = (0:xlen-1)/fs;
%
% wlen = 2^9; % 窗函数宽度取 2^9
% hop = wlen/4;
% nfft = 4*wlen;
% % win = window(@blackman,wlen);
% win = blackmanharris(wlen, 'periodic');
% %% FFT
% figure(NumberTitle="off",Name="音频信号频域图像");
% X = fft(x,nfft);
% f = linspace(-fs/2,fs/2,nfft);
% plot(f,abs(X))
% xlabel("频率f"); ylabel("幅度A"); title('音频信号频域图像')
% %% STFT
% [S1, f1, t1] = stft(x, fs, 'Window', win, 'OverlapLength',wlen - hop,'FFTLength',nfft);
% figure(NumberTitle="off",Name="时频分析图");
% PlotSTFT(t1, f1, S1);
%% example2
% chirp信号仿真
T = 100e-9; B = 77e9; fs = 4*B; c = 3e8;
N = round(T*fs);
t = linspace(0,T,N);
f = linspace(-fs/2,fs/2,N);
k = B/T;
% chirp 信号定义
chirp = exp(1i*pi*k*t.^2);
% 发射信号绘制
% figure(NumberTitle="off",Name='chirp信号绘制')
% subplot(2,1,1)
% plot(t.*1e9,real(chirp));
% xlabel('时间 ns'); ylabel('幅度'); title('chirp 信号 时域')
% subplot(2,1,2)
% plot(f./1e9,fftshift(abs(fft(chirp)/N)))
% xlabel('频率 MHz'); ylabel('幅度'); title('chirp 信号 频域')
% 增加目标
R = [1 5 10];
Sr = zeros(1,N);
delta_t = 2*R/c;
for i = 1:length(R)
Sr = Sr + exp(1i*pi*k*(t - delta_t(i)).^2);
end
Sr = chirp .* conj(Sr);% 匹配滤波
% 接收回波
figure(NumberTitle="off",Name='接收回波')
subplot(3,1,1)
plot(t.*1e9,real(Sr))
xlabel('时间 ns'); ylabel('幅度'); title('接收信号 时域')
dsit = f*c*T/2/B;
subplot(3,1,2)
plot(dsit,fftshift(abs(fft(Sr)/N)))
xlabel('距离 m'); ylabel('幅度'); title('接收信号 距离图')
subplot(3,1,3)
plot(f./1e9,fftshift(abs(fft(Sr)/N)))
xlabel('频率 MKz'); ylabel('幅度'); title('接收信号 频域')
%% STFT
%==========================================================================
% 满足 Constant OverLap-Add STFT
wlen = 2^8;
noverlap = wlen/2; % 步长为窗的长度
% 构造满足 Constant OverLap-Add 长度矩阵
y = zeros(1,wlen - mod(N,noverlap) + N);
y(1:N) = Sr;
% 增加点数后信号参数：
N1 = length(y); T1 = N1/fs;
filter_win = zeros(1,N1); % 初始化窗矩阵
% win = hann(wlen,'periodic');
win = hamming(wlen);
loop = (N1 - wlen)/noverlap;
for i = 0:loop
filter_win(1 + i*noverlap:wlen + i*noverlap) = filter_win(1 + i*noverlap:wlen + i*noverlap) + win';
end
tt = linspace(0,T1,N1);
figure(NumberTitle="off",Name='采样点序列 wlen=N/2 时的COLA验证')
plot(tt.*1e9,filter_win)
xlabel('时间 ns'); ylabel('幅度'); title('采样点序列')
%==========================================================================
% STFT
nfft = length(y);
[S1, f1, t1] = stft(y, fs, 'Window', win, 'OverlapLength', noverlap, 'FFTLength',nfft);
figure(NumberTitle="off",Name="时频分析图");
PlotSTFT(t1, f1, S1);
%% STFT信号重构
x2 = istft(S1,fs,'Window',win,'OverlapLength',noverlap,'FFTLength',nfft);
figure(NumberTitle="off",Name='STFT信号重构 istft函数')
plot(tt.*1e9,real(x2))
xlabel('时间 ns'); ylabel('幅度'); title('重构序列')
