【⾳频处理】之Matlab实现信号的时域和频域的滤波⽂章⽬录
1. 频域滤波
% ccrma.stanford.edu/~jos/parshl/Overlap_Add_Synthesis.html
function overlap_add_synthesis
clo all;
infile = 'P501_C_chine_f1_FB_48k.wav';
infile = 'P501_C_chine_m1_FB_48k.wav';
outfile = regexprep(infile, '\.wav', '_20200527_new.wav');
[data, fs] = audioread(infile); % 384000, 48k
%%
figure;
plot(data);
title('原始信号时域');
xlabel('t/s');
ylabel('幅度');
%%
frame_length = 1024;
frame_shift = 512;
frame_overlap = frame_length - frame_shift;
fftsize = 2^nextpow2(frame_length);
组合图片% hamming windows: sqrt would make the values larger, periodic
analysis_window = sqrt(hamming(frame_length, 'periodic'));
% analysis_window = ones(frame_length, 1); square windows
synthesis_window = analysis_window;
% 1024 * 749 (frame_length * num_frames)
x = buffer(data, frame_length, frame_overlap, 'nodelay'); % 'nodelay' means start from zero, withour padding
% Analysis windowing: (1024, 749) * (1024,1) could be multiplied by .*
x_anawin = x.*repmat(analysis_window(:), [1, size(x, 2)]);
% To frequency-domain
X = fft(x_anawin, fftsize); % default: fft(x_anawin, fftsize, 1);
figure;
一张老照片f = linspace(-24000, 24000, 1024);
plot(abs(X));
title('原始信号频域');
xlabel('f/Hz');
ylabel('幅度');桐油有毒吗
%%
% Modify the spectrum X here
% generate filter, i.e., (1024, 1)
x_f=[0, 1000, 2000, 3500, 24000];
y_f=[1, 0.98, 0.7, 0.55, 0.1];
% 1 means DC offt, X(k) should equal X(1026-k), thus X(513) is the
% center, X(512) = X(514), X(2) = X(1024), X(1) = DC offt.
xx=1:48000/1024:48000/1024*513;
fil=interp1(x_f,y_f,xx,'cubic');
% [513, 1026] 513 points, while X(513) and X(last) should be removed,
% since X(513) is contained in fil, X(last) is not ud(X(0) has no synmetric point)
fil2=fliplr(fil);
fil2(1) = [];
fil2(512) = [];
fil = [fil fil2];
figure;
plot(fil);
X = X.* fil';
% X = X.* 0.2;
figure;
plot(f, abs(X));
% X(1:513,:) = X(1:513,:) * 0.2;
%%
% To time-domain
y = ifft(X);
% Synthesis windowing
y_synwin = y.* repmat(synthesis_window(:), [1, size(y, 2)]);
% Overlap-add
y_ola = y_synwin(:, 1);
for i = 2:size(y_synwin, 2)
y_ola(end+(1:frame_shift)) = zeros(frame_shift, 1);
y_ola(end+(-frame_length+1:0)) = y_ola(end+(-frame_length+1:0)) + y_synwin(:, i); end
% Output
audiowrite(outfile, y_ola, fs);
end
司马穰苴
对应的原始时域信号、频域信号、滤波器、滤波后的频域信号分别为：
⽤adobe audition查看，即可看到频率被减⼩了，即实现了⾼通滤波器。注：红⾊为原来的频率，黄⾊为滤波后的。
2. 时域滤波
f1=5;%第⼀个点频信号分量频率
f2=15;%第⼆个点频信号分量频率
fs=100;%采样率
T=2;%时宽
B=10;%FIR截⽌频率
B=10;%FIR截⽌频率
n=round(T*fs);%采样点个数
t=linspace(0,T,n);
y1= cos(2*pi*f1*t);
y2=cos(2*pi*f2*t);
y=y1+y2;
%%
figure;
plot(t,y);
title('原始信号时域');
xlabel('t/s');
ylabel('幅度');
figure;
plot(t,y1);
title('原始信号时域');
xlabel('t/s');
ylabel('幅度');
figure;
plot(t,y2);
title('原始信号时域');
xlabel('t/s');
ylabel('幅度');
%%
figure;
fft_y=fftshift(fft(y));
f=linspace(-fs/2,fs/2,n);
plot(f,abs(fft_y));
title('原始信号频谱');
xlabel('f/Hz');
ylabel('幅度');
axis([ 0 50 0 100]);家常菜好吃
figure;
fft_y=fftshift(fft(y1));
f=linspace(-fs/2,fs/2,n);果腊
plot(f,abs(fft_y));
title('原始信号y1频谱');
xlabel('f/Hz');
ylabel('幅度');音级
axis([ 0 50 0 100]);
figure;
fft_y=fftshift(fft(y2));
f=linspace(-fs/2,fs/2,n);
plot(f,abs(fft_y));
title('原始信号y2频谱');
xlabel('f/Hz');
ylabel('幅度');
axis([ 0 50 0 100]);
%%
b=fir1(80, B/(fs/2),'high'); %滤波产⽣指定带宽的噪声信号figure;
freqz(b);
%%拨浪鼓图片
y_after_fir=filter(b,1,y);
%%
figure;
plot(t,y_after_fir);
title('滤波后信号时域');
