using DSP

# DSP.jl 提供的主要功能：
# - FFT / IFFT
# - 窗函数（Window Functions）
# - FIR / IIR 滤波器设计
# - 卷积
# - 频谱分析
# - 重采样

# 同时使用 FFTW（DSP.jl 的 FFT 后端）
using FFTW

using FFTW

# 生成信号：两个正弦波叠加
fs = 1000            # 采样率 1000 Hz
t = 0:1/fs:1-1/fs    # 1 秒时间序列
f1, f2 = 50, 120     # 50 Hz 和 120 Hz
signal = sin.(2π * f1 .* t) + 0.5 * sin.(2π * f2 .* t)
signal .+= 0.3 * randn(length(t))    # 添加噪声

# FFT
Y = fft(signal)

# 频率轴
freqs = fftfreq(length(t), fs)    # 自动处理频率排列

# 功率谱
power = abs2.(Y) / length(t)^2

# 只取正频率部分
n = length(t) ÷ 2
freqs_pos = freqs[1:n]
power_pos = 2 * power[1:n]       # ×2 补偿负频率

# IFFT（逆变换）
signal_reconstructed = real(ifft(Y))

using Plots

# 绘制时域信号
p1 = plot(t[1:200], signal[1:200],
    xlabel="时间 (s)", ylabel="振幅",
    title="时域信号", label="原始信号")

# 绘制频谱
p2 = plot(freqs_pos, power_pos,
    xlabel="频率 (Hz)", ylabel="功率",
    title="频谱", label="功率谱",
    xlims=(0, 200))

plot(p1, p2, layout=(2,1), size=(800, 500))

using DSP.Windows

N = 64    # 窗长度

# 常用窗函数
rect_win = rect(N)               # 矩形窗
hann_win = hanning(N)            # Hanning 窗
hamming_win = hamming(N)         # Hamming 窗
blackman_win = blackman(N)       # Blackman 窗
kaiser_win = kaiser(N, 5.0)     # Kaiser 窗（β=5）
gauss_win = gaussian(N, 0.4)    # 高斯窗
tukey_win = tukey(N, 0.5)       # Tukey 窗

using Plots

# 窗函数对比
N = 64
plot(rect(N), label="矩形", linewidth=2)
plot!(hanning(N), label="Hanning", linewidth=2)
plot!(hamming(N), label="Hamming", linewidth=2)
plot!(blackman(N), label="Blackman", linewidth=2)
plot!(kaiser(N, 5.0), label="Kaiser(β=5)", linewidth=2)
xlabel!("样本点"); ylabel!("幅度"); title!("窗函数对比")

using DSP, DSP.Periodograms

fs = 1000
t = 0:1/fs:5-1/fs
signal = sin.(2π * 50 .* t) + 0.5 * sin.(2π * 120 .* t) + randn(length(t)) * 0.3

# 周期图（Periodogram）
pgram = periodogram(signal; fs=fs)
freqs_p = freq(pgram)
power_p = power(pgram)

# Welch 方法（分段平均，减少方差）
pgram_welch = welch_pgram(signal; fs=fs, nfft=256)
freqs_w = freq(pgram_welch)
power_w = power(pgram_welch)

# 短时傅里叶变换（STFT）
S = spectrogram(signal, 256, 128; fs=fs)    # 窗长 256, 重叠 128
t_stft = time(S)
f_stft = freq(S)
S_power = power(S)

using Plots

# Welch 频谱
plot(freqs_w, 10 .* log10.(power_w),
    xlabel="频率 (Hz)", ylabel="功率 (dB)",
    title="Welch 功率谱密度", label="PSD",
    xlims=(0, 200))

# 时频图（Spectrogram）
heatmap(t_stft, f_stft, 10 .* log10.(S_power'),
    xlabel="时间 (s)", ylabel="频率 (Hz)",
    title="频谱图", color=:viridis)

using DSP

# 低通 FIR 滤波器
fs = 1000
f_cutoff = 100    # 截止频率 100 Hz
N = 64            # 滤波器阶数

# 窗函数法设计
lpf = digitalfilter(Lowpass(f_cutoff; fs=fs), FIRWindow(hamming(N)))
# lpf 是滤波器系数向量

# 高通滤波器
hpf = digitalfilter(Highpass(200; fs=fs), FIRWindow(hamming(N)))

# 带通滤波器
bpf = digitalfilter(Bandpass(50, 200; fs=fs), FIRWindow(hamming(N)))

# 带阻滤波器
bsf = digitalfilter(Bandstop(80, 120; fs=fs), FIRWindow(hamming(N)))

# 滤波
filtered_signal = filt(lpf, signal)

# 频率响应
H = freqresp(PolynomialRatio(lpf, [1.0]))

using DSP

fs = 1000

# Butterworth 低通
lpf_butter = digitalfilter(Lowpass(100; fs=fs), Butterworth(4))   # 4 阶

# Chebyshev Type I（通带纹波）
lpf_cheby1 = digitalfilter(Lowpass(100; fs=fs), Chebyshev1(4, 1))   # 4 阶, 1dB 纹波

# Chebyshev Type II（阻带纹波）
lpf_cheby2 = digitalfilter(Lowpass(100; fs=fs), Chebyshev2(4, 40))  # 4 阶, 40dB 衰减

# 椭圆滤波器
lpf_ellip = digitalfilter(Lowpass(100; fs=fs), Elliptic(4, 1, 40))  # 4 阶, 1dB/40dB

# 应用滤波
y = filt(lpf_butter, signal)

# 二阶节级联（更稳定）
sos = SecondOrderSections(lpf_butter)
y_stable = filt(sos, signal)

using DSP

# 线性卷积
x = [1, 2, 3, 4]
h = [1, 0, -1]
y = conv(x, h)      # [1, 2, 2, 2, -4, -4]

# 2D 卷积（图像处理）
img = rand(100, 100)
kernel = ones(3, 3) / 9    # 3×3 均值滤波
filtered = conv(img, kernel)

# 互相关（用于延迟检测）
x = randn(100)
y = circshift(x, 5) + randn(100) * 0.1    # 延迟 5 个样本
r = xcorr(y, x)
delay = argmax(r) - length(x)

using DSP

# 重采样（改变采样率）
fs_old = 44100       # CD 品质
fs_new = 16000       # 语音处理常用

# 上采样 + 抗混叠滤波 + 下采样
ratio = gcd(fs_old, fs_new)
up = fs_new ÷ ratio
down = fs_old ÷ ratio

signal_44k = sin.(2π * 440 .* (0:1/fs_old:1-1/fs_old))
signal_16k = resample(signal_44k, up // down)

# 插值
using Interpolations

t_orig = 0:0.1:5
y_orig = sin.(t_orig)

# 线性插值
itp_linear = LinearInterpolation(t_orig, y_orig)
y_fine = itp_linear.(0:0.01:5)

# 三次样条插值
itp_cubic = CubicSplineInterpolation(t_orig, y_orig)
y_fine_cubic = itp_cubic.(0:0.01:5)

using DSP

# STFT 参数
fs = 1000
t = 0:1/fs:2-1/fs

# 生成变频信号（chirp）
f_inst = 50 .+ 200 .* t     # 瞬时频率从 50 Hz 到 250 Hz
signal = sin.(2π .* cumsum(f_inst) / fs)

# STFT
window = hamming(256)
S = spectrogram(signal, 256, 192; fs=fs, window=window)

# 绘制时频图
using Plots
heatmap(time(S), freq(S), 10 .* log10.(power(S)'),
    xlabel="时间 (s)", ylabel="频率 (Hz)",
    title="Chirp 信号 STFT",
    color=:inferno)

# 从 STFT 重建信号
function istft(S, hop; nfft=size(S, 1))
    n_frames = size(S, 2)
    signal_len = (n_frames - 1) * hop + nfft
    signal = zeros(signal_len)
    window_sum = zeros(signal_len)
    w = hamming(nfft)
    for i in 1:n_frames
        idx = (i - 1) * hop + 1
        signal[idx:idx+nfft-1] .+= real(ifft(S[:, i]))[1:nfft] .* w
        window_sum[idx:idx+nfft-1] .+= w .^ 2
    end
    return signal ./ max.(window_sum, 1e-10)
end

using WAV, DSP, FFTW

# 读取音频文件
audio, fs = wavread("audio.wav")

# 如果是立体声，取单声道
if ndims(audio) > 1
    audio = audio[:, 1]
end

# 降噪：频域门限法
Y = fft(audio)
power = abs2.(Y)

# 设定门限
threshold = median(power) * 3
Y_denoised = copy(Y)
Y_denoised[power .< threshold] .= 0    # 抑制低功率成分

audio_clean = real(ifft(Y_denoised))

# 低通滤波去除高频噪声
lpf = digitalfilter(Lowpass(4000; fs=fs), Butterworth(6))
audio_filtered = filt(lpf, audio_clean)

# 保存
wavwrite(audio_filtered, "audio_clean.wav"; Fs=fs)

using DSP, FFTW

# AM 调幅信号解调
fs = 10000
t = 0:1/fs:0.1-1/fs

# 消息信号
msg = sin.(2π * 100 .* t)

# 载波信号
carrier_freq = 1000
carrier = cos.(2π * carrier_freq .* t)

# AM 调制
m = 0.5    # 调制指数
am_signal = (1 .+ m .* msg) .* carrier

# 包络检波（简单解调）
# 全波整流 + 低通滤波
rectified = abs.(am_signal)
lpf = digitalfilter(Lowpass(200; fs=fs), Butterworth(4))
demodulated = filt(lpf, rectified)

# 绘制
using Plots
p1 = plot(t[1:500], msg[1:500], title="消息信号", label="msg")
p2 = plot(t[1:500], am_signal[1:500], title="AM 信号", label="AM")
p3 = plot(t[1:500], demodulated[1:500], title="解调信号", label="解调")
plot(p1, p2, p3, layout=(3,1), size=(800, 600))

using DSP, FFTW

# 模拟轴承振动信号
fs = 10000
t = 0:1/fs:2-1/fs

# 正常信号 + 故障特征频率
normal = 0.1 * randn(length(t))
fault_freq = 150    # 故障特征频率 150 Hz
fault = 0.5 * sin.(2π * fault_freq .* t) .* (mod.(t, 0.01) .< 0.001)

signal = normal + fault

# 频谱分析
pgram = welch_pgram(signal; fs=fs, nfft=1024)

# 包络分析（用于轴承故障检测）
analytic = hilbert(signal)         # 解析信号
envelope = abs.(analytic)          # 包络
env_spectrum = welch_pgram(envelope; fs=fs, nfft=1024)

using Plots
p1 = plot(freq(pgram), 10 .* log10.(power(pgram)),
    title="频谱", label="PSD", xlims=(0, 500))
p2 = plot(freq(env_spectrum), 10 .* log10.(power(env_spectrum)),
    title="包络谱", label="包络 PSD", xlims=(0, 500))
plot(p1, p2, layout=(2,1), size=(800, 500))

using Wavelets

# 离散小波变换（DWT）
signal = randn(128)

# Daubechies 小波
wt = wavelet(WT.db4)    # 4 阶 Daubechies
coeffs = dwt(signal, wt)

# 逆变换
signal_recon = idwt(coeffs, wt)

# 多层分解
level = 3
coeffs_multi = dwt(signal, wt, level)

# 小波去噪
# 软阈值法
function wavelet_denoise(signal, wt; level=4, threshold_factor=1.5)
    coeffs = dwt(signal, wt, level)
    # 估计噪声标准差
    σ = median(abs.(coeffs)) / 0.6745
    threshold = threshold_factor * σ
    # 软阈值
    coeffs_thresh = sign.(coeffs) .* max.(abs.(coeffs) .- threshold, 0)
    return idwt(coeffs_thresh, wt, level)
end