请用matlab代码利用谱减法消除一段带有噪声的信号

假设你有一个带噪声的音频信号，可以按照以下步骤使用谱减法进行降噪。

1. 读取音频信号

[x, fs] = audioread('noisy_audio.wav');

2. 计算短时傅里叶变换（STFT）并获取幅度谱和相位谱

win_length = round(fs * 0.02); % 窗口长度为20ms
hop_length = round(fs * 0.01); % 窗口重叠为10ms
n_fft = 2^nextpow2(win_length); % FFT点数
[S, f, t] = spectrogram(x, hamming(win_length), hop_length, n_fft, fs); % STFT
A = abs(S); % 幅度谱
P = angle(S); % 相位谱

3. 计算噪声的幅度谱

n_fft2 = 2^nextpow2(win_length / 2 + 1); % 计算FFT点数（取一半的长度加1）
N = zeros(n_fft2, 1); % 初始化噪声幅度谱
for i = 1:size(A, 2)
    N = N + A(:, i).^2; % 累加所有帧的幅度谱的平方
end
N = N / size(A, 2); % 求平均值

4. 估计信号的幅度谱

alpha = 2; % 设定谱减系数
A_est = max(A - alpha * sqrt(N), 0); % 谱减

5. 合成信号

S_est = A_est .* exp(1j * P); % 合成幅度谱和相位谱
x_est = istft(S_est, hamming(win_length), hop_length, n_fft, fs); % 合成信号

其中，`istft` 是用于反向转换的函数，可以使用以下代码定义：

function x = istft(S, win, hop, n_fft, fs)
% ISTFT Inverse short-time Fourier transform
%   x = ISTFT(S, WIN, HOP, N_FFT, FS) computes the inverse short-time
%   Fourier transform (ISTFT) of the STFT coefficients S using a window
%   function WIN with hop size HOP and FFT length N_FFT. The sampling rate
%   is FS. The output signal x has the same length as the original signal.
%
%   Inputs:
%       S:      N_FFT/2+1 x T matrix of STFT coefficients
%       WIN:    N x 1 window function (N <= N_FFT)
%       HOP:    hop size in samples
%       N_FFT:  FFT length
%       FS:     sampling rate in Hz
%
%   Outputs:
%       x:      L x 1 vector of reconstructed signal

% Get the number of frequency bins and frames
[n_bins, n_frames] = size(S);

% Check the window size
win_len = length(win);
if win_len > n_fft
    error('Window size exceeds FFT length.');
end

% Initialize the reconstructed signal
L = (n_frames - 1) * hop + win_len;
x = zeros(L, 1);

% Compute the synthesis window
win_syn = win / sum(win);

% Synthesize each frame
for i = 1:n_frames
    % Get the current frame
    s = S(:, i);
    
    % Compute the magnitude spectrum
    a = [s; conj(s(end-1:-1:2))]; % reflect upper half of the spectrum
    x_mag = real(ifft(a));
    x_mag = x_mag(1:win_len);
    
    % Overlap-add the reconstructed frame
    idx = (i-1)*hop + (1:win_len);
    x(idx) = x(idx) + x_mag .* win_syn;
end

% Normalize the reconstructed signal
x = x / max(abs(x));
end