数字信号处理

Posted by Jason Zhu Blog on January 30, 2024

1.计算正弦信号的频谱

close all; clear all;
addpath 'pdm-0925'
addpath 'U:\cvr\pdm\pdm_fix\source\'


file_path = 'E:\matlab\pdm\pdmmodule\pdm_new-1011\pdm_16x48k_test5\';

fid = fopen([file_path 'fifo_16bit.pcm'],'rb');
Sig = fread(fid,'bit16');
fclose(fid);
data = double(Sig) / 2^16;
fs = 48000;
sig_fs = 997;

test_start = 0.8 * fs;
test_end =  1.8*fs -1;

snr(data(test_start : test_end, 1), fs);
figure(2);

ft = fft(data(test_start : test_end, 1));
L = test_end - test_start;
P2 = abs(ft/L);
P1 = P2(1:L/2+1);%单边谱

plot( 20 * log10(abs(P1)));
%audiowrite('swap_32000_16bit.wav', data, fs);
% [Sig, fs] = audioread('pdm-16kfs-1k-16bit-3.072MHz-0925.wav');
% sinad(Sig(2*fs:3*fs-1, 1), fs);

2.数据截位仿真

clc;
a = 10000 *randn(1,1024);
fs = -512 : 1 : 511;
b = floor(a / 2^5);
afft = fftshift(fft(a));
bfft = fftshift(fft(b));

figure(1);
subplot(2, 1, 1);
plot(fs, abs(afft));
title('原始信号');
subplot(2, 1, 2);
plot(fs, abs(bfft));
title('截位信号');
b = floor(a/2^12); %截去低12位
afft = fftshift(fft(a));
bfft = fftshift(fft(b));
figure(2);
subplot (2,1,1);
plot(fs, abs(afft));
title('原始信号');

subplot (2,1,2);
plot(fs, abs(bfft));
title('截位信号');


b = zeros(1,length(a)); %截去低12位
for ii = 1:length(a)
    if a(ii) < 0
        b(ii) = floor(a(ii)/2^12) + 1;
    else
        b(ii) = floor(a(ii)/2^12);
end
end

afft = fftshift(fft(a));
bfft = fftshift(fft(b));
figure(3);
subplot(2,1,1);
plot(fs, abs(afft));
title('原始信号');
subplot(2, 1, 2);
plot(fs, abs(bfft));
title('截位信号');



b = zeros(1,length(a)); %截去低12位
for ii = 1:length(a)
    if a(ii) < 0
        b(ii) = -floor(-a(ii)/2^12);% 负数取绝对值后再截低12位再变回负数
    else
        b(ii) = floor(a(ii)/2^12); % 直接截位
    end
end
afft = fftshift(fft(a) );
bfft = fftshift(fft(b));
figure(4);
subplot(2,1,1);
plot(fs, abs(afft));
title('原始信号');
subplot(2, 1, 2);
plot(fs, abs(bfft));
title('截位信号');

1bit截取误差：

clc;
a = 10000 *randn(1,10240);
b = floor(a / 2^12) / 32768;

figure(3);
L=10240;
fs_l = 0 : 1 : 10239;
ft = fft(b);
P2 = abs(ft/L);
P1 = P2(1:L/2+1);
plot( 20 * log10(abs(P1 * 2)));





clc;
clear all;
close all;

f = 997; % 1 kHz的频率
fs = 48000; % 48 kHz的采样率
t = 0:1/fs:1; % 从0到1秒的时间矢量
a = (sin(2 * pi * f * t) - 0.5) * 2^23;
x = bitand(int32(a), 255) - 128;
b = double(floor(double(x) / 2^8)) / 32768;

figure(3);
L=48001;
fs_l = 0 : 1 : 48000;
ft = fft(b);
P2 = abs(ft/L);
P1 = P2(1:L/2+1);
plot( 20 * log10(abs(P1 * 2)));

3.自适应滤波仿真

% 生成一个原始信号
Fs = 1000; % 采样频率
t = 0:1/Fs:1; % 时间向量
f = 5; % 信号频率
original_signal = sin(2*pi*f*t);

% 添加噪声
N = 9;
square_wave = zeros(size(t));
for n = 3:2:N
    square_wave = square_wave + (4/pi) * (1/n) * sin(2*pi*n*f*t);
end
noisy_signal = original_signal + square_wave;

% 创建自适应滤波器
filter_order = 10; % 滤波器阶数
mu = 0.01; % 步长（学习率）
adaptive_filter = dsp.LMSFilter('StepSize', mu, 'Length', filter_order);

% 初始化输出和误差存储
filtered_signal = zeros(size(noisy_signal));
error_signal = zeros(size(noisy_signal));

% 自适应滤波
for i = 1:length(noisy_signal)
    [filtered_signal(i), error_signal(i)] = step(adaptive_filter, noisy_signal(i), original_signal(i));
end

% 绘制原始信号、带噪声信号和滤波后信号
figure;
subplot(3,1,1);
plot(t, original_signal);
title('原始信号');

subplot(3,1,2);
plot(t, noisy_signal);
title('带噪声信号');

subplot(3,1,3);
plot(t, filtered_signal);
title('滤波后信号');

% 绘制误差信号
figure;
plot(t, error_signal);
title('误差信号');



%%%%%%%%%%%%%%FFT%%%%%%%%%%%%%%%
ft_o = fft(noisy_signal);
L = 1001;
P2_o = abs(ft_o/L);
P1_o = P2_o(1:L/2+1);
figure;
plot( 20 * log10(abs(P1_o)));

ft_o = fft(filtered_signal);
L = 1001;
P2_o = abs(ft_o/L);
P1_o = P2_o(1:L/2+1);
figure;
plot( 20 * log10(abs(P1_o)));

自适应滤波仿真

clc;%清除工作区
close all;
clear all;%清除命令窗口

t = 0.1:0.1:100; %时间轴
signal_len = length(t);%信号长度
d = sin(t);%理想信号

noise = 0.1*randn(1,signal_len);%噪声
d_noise =  d+noise;%要过滤的信号-加了噪声的信号

filter_len = 50;  %滤波器长度
W = zeros(1,filter_len); %初始化滤波器
x = zeros(1,filter_len); %初始化卷积输入
After_filter = zeros(1,signal_len);%用于存储滤波之后的信号
e = zeros(1,signal_len); %初始化误差
mull = 0.03; %步长


for k =1:signal_len
    x = [d_noise(k) x(1:filter_len-1)]; %线性卷积的输入
    After_filter(k) = W*x';
    e(k) =d(k) - After_filter(k);%计算误差
    W = W + 2*mull*e(k)*x; %计算误差更新权重
end

figure(1)
subplot(411)
plot(t,d)
title('LMS')
xlabel('理想信号')

subplot(412)
plot(t,d_noise)
xlabel('含噪原信号')

subplot(413)
plot(t,e)
xlabel('误差')

subplot(414)
plot(t,After_filter)
xlabel('去噪信号')

自适应滤波仿真

clc;%清除工作区
close all;
clear all;%清除命令窗口

Fs = 48000;
t = 0 :1 / Fs:1; %时间轴
signal_len = length(t);%信号长度
f = 997; 
d = sin(2*pi*f*t);%理想信号

noise = 0.01*randn(1,signal_len);%噪声
d_noise =  d + noise;%要过滤的信号-加了噪声的信号

filter_len = 50;  %滤波器长度
W = zeros(1,filter_len); %初始化滤波器
x = zeros(1,filter_len); %初始化卷积输入
After_filter = zeros(1,signal_len);%用于存储滤波之后的信号
e = zeros(1,signal_len); %初始化误差
mull = 0.03; %步长


for k =1:signal_len
    x = [d_noise(k) x(1:filter_len-1)]; %线性卷积的输入
    After_filter(k) = W*x';
    e(k) =d(k) - After_filter(k);%计算误差
    W = W + 2*mull*e(k)*x; %计算误差更新权重
end

figure(1)
subplot(411)
plot(t,d)
title('LMS')
xlabel('理想信号')

subplot(412)
plot(t,d_noise)
xlabel('含噪原信号')

subplot(413)
plot(t,e)
xlabel('误差')

subplot(414)
plot(t,After_filter)
xlabel('去噪信号')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%频谱%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ft_o = fft(d);
L = signal_len;
P2_o = abs(ft_o/L);
P1_o = P2_o(1:L/2+1);
figure(2);
plot( 20 * log10(abs(2* P1_o)));
title('d fft');


ft_o = fft(d_noise);
L = signal_len;
P2_o = abs(ft_o/L);
P1_o = P2_o(1:L/2+1);
figure(3);
plot( 20 * log10(abs(2* P1_o)));
title('d_noise fft');


ft_o = fft(After_filter);
L = signal_len;
P2_o = abs(ft_o/L);
P1_o = P2_o(1:L/2+1);
figure(4);
plot( 20 * log10(abs(2* P1_o)));
title('After_filter fft');


figure(5);
snr(d,Fs);
title('snr d fft');

figure(6);
snr(d_noise,Fs);
title('snr d_noise fft');

figure(7);
snr(After_filter,Fs);
title('snr After_filter fft');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%signal 2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f_2 = 4997; 
d_2 = sin(2*pi*f_2*t);%理想信号

N = 9;
square_wave_2 = zeros(size(t));

for n = 3:2:N
    square_wave_2 = square_wave_2 + (4/pi/100) * (1/n) * sin(2*pi*n*f*t);
end
d_noise_2 =  d_2 + noise + square_wave_2;%要过滤的信号-加了噪声的信号


filtered_signal = filter(W, 1, d_noise_2);
figure(8);
snr(d_noise_2,Fs);
title('snr d_noise_2 fft');
figure(9);
snr(filtered_signal,Fs);
title('snr filtered_signal fft');

自适应滤波

% 生成原始正弦信号
clc;
clear all;
close all;
Fs = 1000;               % 采样频率
t = 0:1/Fs:1-1/Fs;       % 时间向量
f = 50;                  % 正弦信号频率
A = 1;                   % 正弦信号幅度
cleanSignal = A * sin(2*pi*f*t);  % 原始正弦信号

% 生成已知的噪声信号
noise = 0.001 * randn(size(t));     % 添加高斯噪声

% 带噪声信号
noisySignal = cleanSignal + noise;

% 创建自适应滤波器对象
mu = 0.01;               % 学习速率
filterLength = 32;      % 滤波器长度
adaptFilter = dsp.LMSFilter('StepSize', mu, 'Length', filterLength);

% 初始化滤波后的信号
filteredSignal = zeros(size(noisySignal));

% 使用自适应滤波器来滤除噪声
for i = 1:length(noisySignal)
    [y, ~] = adaptFilter(noisySignal(i), noise(i));  % 使用已知噪声信号作为参考信号
    filteredSignal(i) = y;
end

% 绘制结果，比较原始信号、带噪声信号和滤波后的信号
figure;
subplot(3,1,1);
plot(t, cleanSignal);
title('原始正弦信号');

subplot(3,1,2);
plot(t, noisySignal);
title('带噪声正弦信号');

subplot(3,1,3);
plot(t, filteredSignal);
title('滤波后的正弦信号');

xlabel('时间 (秒)');

filter_s = noisySignal - filteredSignal;
figure;
plot(t, filter_s);

figure;
snr(noisySignal, Fs);

figure;
snr(filter_s, Fs);

4.卡尔曼滤波

clc;
clear all;
close all;

[ys, Fs] = audioread("pdm-48k-16bit-0db.wav");

noisy_signal = ys(40000:50000,1);

% 初始化卡尔曼滤波器参数
A = 1;  % 状态转移矩阵
H = 1;  % 观测矩阵
Q = 0.001;  % 系统过程噪声协方差
R = 0.1;  % 观测噪声协方差
x = 0.5;  % 初始振幅估计
P = 1;  % 初始估计误差协方差

% 初始化估计的振幅数组
estimated_amplitudes = zeros(size(noisy_signal));

% 卡尔曼滤波过程
for i = 1:length(noisy_signal)
    % 预测步骤（振幅估计）
    x_hat = A * x;
    P = A * P * A' + Q;
    
    % 更新步骤（振幅校正）
    K = P * H' / (H * P * H' + R);
    x = x_hat + K * (noisy_signal(i) - H * x_hat);
    P = (1 - K * H) * P;
    
    % 保存估计的振幅
    estimated_amplitudes(i) = x;
end

% 绘制原始信号、噪声信号和估计的振幅
figure;
subplot(3,1,2);
plot(noisy_signal, 'g');
title('带噪声的音频信号');
xlabel('时间 (秒)');
subplot(3,1,3);
plot(estimated_amplitudes, 'r');
title('估计的振幅');
xlabel('时间 (秒)');

figure;
snr(noisy_signal,Fs);
figure;
snr(estimated_amplitudes,Fs);

5.自动增益控制

clc;
clear all;
close all;
% 1. 读取音频文件
[inputAudio, Fs] = audioread('rk3588-48k_96db.wav');

% 2. 设置AGC参数
targetRMS = 0.5; % 目标RMS值
attackTime = 0.01; % 增益调整的攻击时间（秒）
releaseTime = 0.1; % 增益调整的释放时间（秒）

% 3. 初始化增益和平滑增益
gain = ones(size(inputAudio));
smoothedGain = ones(size(inputAudio));

% 4. 处理音频数据
for i = 1:length(inputAudio)
    % 4.1 计算瞬时RMS值
    windowSize = round(Fs * 0.02); % 20毫秒窗口大小
    startIdx = max(1, i - windowSize + 1);
    instantRMS = rms(inputAudio(startIdx:i));
    
    % 4.2 计算误差
    error = targetRMS - instantRMS;
    
    % 4.3 计算增益
    if error > 0
        alpha = exp(-1 / (Fs * attackTime));
        gain(i) = gain(i) * (1 - alpha) + alpha;
    else
        beta = exp(-1 / (Fs * releaseTime));
        gain(i) = gain(i) * (1 - beta) + beta;
    end
    
    % 4.4 平滑增益处理
    smoothedGain(i) = max(gain(i), smoothedGain(i));
    
    % 4.5 将增益应用到音频信号
    outputAudio(i) = inputAudio(i) * smoothedGain(i);
end

% 5. 播放输出音频
sound(outputAudio, Fs);

% 6. 保存输出音频到文件
audiowrite('output_audio_file.wav', outputAudio, Fs);

6.fft计算

clc;
clear all;
close all;
% 生成正弦信号
Fs = 48000; % 采样率为2MHz
T = 1/Fs; % 采样间隔
t = 0:T:1-T; % 时间向量，假设信号长度为1秒
f_signal = 1e3; % 正弦信号的频率为100kHz
amplitude = 1; % 信号幅度

sinusoidal_signal = amplitude * sin(2*pi*f_signal*t);

% 进行FFT
N = length(sinusoidal_signal);
Y = fft(sinusoidal_signal);
P2 = abs(Y/N); % 双侧频谱
P1 = P2(1:N/2+1); % 单侧频谱，保留正频率部分
f = Fs*(0:N/2)/N; % 频率向量

% 计算归一化频谱的dB表示
P1_dB = 10 * log10(P1 * 2);

% 绘制频谱图
figure;
plot(f, P1_dB);
xlabel('频率 (Hz)');
ylabel('归一化幅值 (dB)');
title('正弦信号的归一化频谱（dB）');
FEATURED TAGS
FRIENDS