[y,Fs]=audioread('1.wav') Nsamps = length(y); t = (1/Fs)*(1:Nsamps) %Prepare time data for plot %Do Fourier Transform y_fft = abs(fft(y)); %Retain Magnitude y_fft = y_fft(1:Nsamps/2); %Discard Half of Points f = Fs*(0:Nsamps/2-1)/Nsamps; %Prepare freq data for plot %Plot Sound File in Time Domain figure plot(t, y) xlabel('Time (s)') ylabel('Amplitude') title('语音信号时域波形图') %Plot Sound File in Frequency Domain figure plot(f, y_fft) xlim([0 1000]) xlabel('Frequency (Hz)') ylabel('Amplitude') title('语音信号频谱图') % 设计 300Hz～3400Hz 的带 % 计算声音的时间长度 % 仿真系统采样率 % 仿真系统采样时间点 [wav,fs]=audioread('GDGvoice8000.wav'); t_end=1/fs *length(wav); Fs=50000; t=1/Fs:1/Fs:t_end; % 利用插值函数将音频信号的采样率提升为 Fs=50KHz wav=interp1([1/fs:1/fs:t_end],wav,t,'spline'); 通预滤波器 H(z) [fenzi,fenmu]=butter(6,[300 3400]/(Fs/2)); nt = wgn(1,length(t),0.1); % 噪声 nt=nt/(max(abs(nt))); %归一化噪声 wav_noise = wav + nt; % 对音频信号进行滤波 wav_after = filter(fenzi,fenmu,wav_noise); figure(1); subplot(2,1,1); plot(wav(53550:53750)); title('语音信号时域波形'); axis([0 200 -0.3 0.3]); subplot(2,1,2); psd(wav, 10000, Fs); title('语音信号功率谱密度'); axis([0 25000 -20 10]); figure(2); subplot(2,1,1); plot(wav_noise(53550:53750)); 

title('加噪声后的语音信号时域波形'); axis([0 200 -0.3 0.3]); subplot(2,1,2); psd(wav_noise, 10000, Fs); title('加噪声后的语音信号功率谱密度'); axis([0 25000 -20 10]); figure(3); subplot(2,1,1); plot(wav_after(53550:53750)); title('滤波后的语音信号时域波形'); axis([0 200 -0.3 0.3]); subplot(2,1,2); psd(wav_after, 10000, Fs); title('滤波后的语音信号功率谱密度'); axis([0 25000 -20 10]); sound(wav/max(wav), Fs); sound(wav_after/max(wav_after), Fs); wavwrite(demod_out,Fs,'SSBDemod_OUT.wav'); % 播放解调音频 % 播放解调音频 % 保存输出信号 %对 n 点进行傅里叶变换到频域 [y,fs,bits]=wavread('C:\Documents and Settings\Administrator\桌面\cuocuo.wav'); sound(y,fs) % 回放语音信号 n=length(y) %选取变换的点数 y_p=fft(y,n); f=fs*(0:n/2-1)/n; % 对应点的频率 figure(1) subplot(2,1,1); plot(y); title('原始语音信号采样后时域波形'); xlabel('时间轴') ylabel('幅值 A') subplot(2,1,2); plot(f,abs(y_p(1:n/2))); %语音信号的频谱图 title('原始语音信号采样后频谱图'); %语音信号的时域波形图 fp=1500;fc=1700;As=100;Ap=1; (以上为低通滤波器的性能指标） wc=2*pi*fc/fs; wp=2*pi*fp/fs; wdel=wc-wp; beta=0.112*(As-8.7); N=ceil((As-8)/2.285/wdel); wn= kaiser(N+1,beta); ws=(wp+wc)/2/pi; b=fir1(N,ws,wn); figure(3); freqz(b,1); 

（此前为低通滤波器设计阶段）——接下来为去除噪声信号的程序—— x=fftfilt(b,y_z); X=fft(x,n); figure(4); subplot(2,2,1); plot(f,abs(y_zp)); title('滤波前信号的频谱'); subplot(2,2,2); plot(f,abs(X)); title('滤波后信号频谱'); subplot(2,2,3); plot(y_z); title('滤波前信号的波形') subplot(2,2,4);plot(x); title('滤波后信号的波形') %sound(x,fs,bits) %回放滤波后的音频 clc; [x,FS]=audioread(‘F:\1.wav’); x=x(:,1); N=length(x); n=0:N-1; t=n/FS; figure(1); subplot(211); plot(x); %sound(x,FS,bits); title('语音信号时域波形图') y=fft(x,N); f=(FS/N)*(1:N); subplot(212); plot(f(1:105000),abs(y(1:105000))); title('语音信号频谱图'); fs=1200; N=1000; n=1:N; A=importdata('F:\123.mat'); 

B=[]; for i=401:1400; b=A(i); B=[B b]; end [c,l]=wavedec(B,5,'db4'); a5=appcoef(c,l,'db4',5); d5=detcoef(c,l,5); d4=detcoef(c,l,4); d3=detcoef(c,l,3); d2=detcoef(c,l,2); d1=detcoef(c,l,1); figure(1); subplot(6,1,1);plot(a5); subplot(6,1,2);plot(d5); subplot(6,1,3);plot(d4); subplot(6,1,4);plot(d3); subplot(6,1,5);plot(d2); subplot(6,1,6);plot(d1); thr1=thselect(d5,'rigrsure')%stein 无偏估计阈值; thr1 =0.1688 max(d5) ans =0.1688 thr1=thselect(d4,'rigrsure')%stein 无偏估计阈值; thr1 =0.4431 thr1=thselect(d4,'rigrsure')%stein 无偏估计阈值; thr1 =0.4431 thr1=thselect(d3,'rigrsure')%stein 无偏估计阈值; thr1 =1.0841 thr1=thselect(d2,'rigrsure')%stein 无偏估计阈值; thr1 =1.5006 

thr1=thselect(d1,'rigrsure')%stein 无偏估计阈值; thr1 =1.2141 cd5=wthresh(d5,'s',0.15); figure(2) subplot(6,1,1);plot(cd5); cd4=wthresh(d4,'s',0.3); subplot(6,1,2);plot(cd4); cd3=wthresh(d3,'s',0.3); subplot(6,1,3);plot(cd3); cd2=wthresh(d2,'s',0.5); subplot(6,1,4);plot(cd2); cd1=wthresh(d1,'s',0.8); subplot(6,1,5);plot(cd1); c1=[a5,cd5,cd4,cd3,cd2,cd1]; ss=waverec(c1,l,'db4'); figure(3); subplot(2,1,1); plot(n*1000/fs,ss); xlabel('时间/ms'); ylabel('幅值/（m/s2）'); subplot(2,1,2);plot(linspace(0,fs/2,N),abs(fftshift(fft(ss)))); plot(linspace(fs/2,fs/2,N),abs(fftshift(fft(ss)))); Y = fft(ss,512); Pyy = Y.* conj(Y) / 512; f = 1000*(0:256)/512; plot(f,Pyy(1:257)) ylabel('幅值/（m/s2）') xlabel('频率/Hz')