Đề tài Lọc tín hiệu (Có code matlab đi kèm)

Contents Lời nói đầu Trong cuộc sống hằng ngày… Lý thuyết chung về lọc nhiễu bằng thuật toán Wiener Filter Thuật toán Wiener filter Giới thiệu Wiener filter là thuật toán được sử dụng rộng rãi trong nâng cao chất lượng tiếng nói. Ý tưởng cơ bản của thuật toán Wiener filter là tạo ra tín hiệu tiếng nói sạch bằng cách nén nhiễu. Ước lượng được thực hiện bằng cách hạ thấp sai số bình phương trung bình (Mean Square Error) giữa tín hiệu mong muốn và tín hiệu ước lượng. Nguyên lý cơ bản của Wiener Filtering Giả thiết rằng y[n] là tín hiệu vào đã bị nhiễu, nó là tổng của tín hiệu sạch s[n] và tín hiệu nhiễu n[n]: y[n]=s[n]+n[n] (3.13) Thực hiện biến đổi Fourier rời rạc cả 2 vế,ta được Y[ω]=S[ω]+N[ω] (3.14) Chúng ta có thể biểu diễn Y() dưới dạng phức như sau: Y[ω]= (3.15) Khi đó Y[ω] là biên độ phổ, và ) là pha của tín hiệu đã bị nhiễu. Phổ của tín hiệu nhiễu N(ω) có thể được biểu diễn dạng biên độ và pha: D(ω) = (3.16) Biên độ phổ của nhiễu D(ω) không xác định được, nhưng có thể thay thế bằng giá trị trung bình của nó được tính trong khi không có tiếng nói(tiếng nói bị dừng), và pha của tín hiệu nhiễu có thể thay thế bằng pha của tín hiệu bị nhiễu , việc làm này không ảnh hưởng đến tính dễ nghe của tiếng nói [27], có thể ảnh hưởng đến chất lượng của tiếng nói là làm thay đổi pha của tiếng nói nhưng cũng chỉ vài độ. Ta có thể ước lượng được biên độ của phổ tín hiệu sạch Ŝ(ω) từ Y() bằng một hàm phi tuyến được xác định như sau : G(ω) = Ŝ(ω)/Y(ω) (3.17) G(ω) Với : G(ω) = (3.18) Đặt Priori SNR và Posteriori SNR như sau : (3.19) (3.20) Một khó khăn trong các thuật toán nâng cao chất lượng tiếng nói là ta không có tín hiệu trước tín hiệu sạch s[n] nên ta không thể biết phổ của nó. Do đó ta không thể tính được mà trong các thuật toán nâng cao chất lượng tiếng nói thì là tham số rất cần thiết để ước lượng tín hiệu sạch.Trong các hệ thống nâng cao chất lượng giọng nói có thể ước lượng được và bằng cách cho các thông số thích hợp vào các phương trình sau : (3.21) (3.22) (3.23) Trong đó P[.] là hàm chỉnh lưu bán sóng có dạng như sau: (3.24) Và và chỉ số để tín hiệu tại khoảng thời gian đang xử lý. Trong phương trình nếu cho hệ số ta có thể ước lượng được bằng . Trong thực tế hệ số =0.98 rất tốt cho các tín hiệu có SNR<4dB. Từ phương trình (3.18) và (3.19) có theo Wiener filter như sau: (3.25) Sơ đồ khối của thuật toán Wiener Filtering: Ước lượng ,cập nhật nhiễu Priori SNR DFT Pha của tín hiệu IDF Hàm xử lý giảm nhiễu Y(ω) Tín hiệu bị nhiễu SNRpri Tín hiệu sau khi tăng cường Hình 1 : Sơ đồ khối của thuật toán Wiener Filter Overlap và Adding trong xử lý tiếng nói Phân tích tín hiệu theo từng frame Do tín hiệu cần xử lý là tín hiệu liên tục, nên khi chúng ta biến đổi DFT trực tiếp tín hiệu từ miền thời gian mà không thông qua một quá trình tiền xử lý nào trước đó thì tín hiệu sau khi được biến đổi DFT sẽ biến đổi nhanh, lúc đó chúng ta không thể thực hiện được các thuật toán xử lý triệt nhiễu trong tín hiệu vì khi đó tín hiệu được xem là động. Chính vì vậy, tín hiệu của chúng ta cần phải được phân tích thành những khung tín hiệu(frame) liên tục trong miền thời gian trước khi chuyển sang miền tần số bằng biến đổi FFT. Khi tín hiệu được phân tích thành các frame liên tục, thì trong từng frame, tín hiệu của chúng ta sẽ biến đổi chậm và nó được xem là tĩnh. Nếu tín hiệu được phân tích theo từng frame thì khi đó các thuật toán xử lý triệt nhiễu trong tín hiệu mới có thể thực hiện được một cách hiệu quả. Và cách phân tích tín hiệu của chúng ta là “frame by frame”. Để thực hiện việc phân tích tín hiệu thành các frame, cần sử dụng các loại cửa sổ thích hợp. Ở đây, chúng ta sử dụng cửa sổ Hamming, với N = 256 mẫu trong từng frame : Overlap và Adding Sau khi phân tích tín hiệu thành các frame liên tục trong miền thời gian bằng cửa sổ Hamming, nếu các frame này liên tục với nhau và không theo một điều kiện nào cả thì khi thực hiện biến đổi FFT thì vô tình chúng ta đã làm suy giảm tín hiệu do Hamming là cửa sổ phi tuyến. Nên khi thực hiện phân tích tín hiệu thành các frame thì yêu cầu đặt ra là các frame phải sắp xếp chồng lên nhau, gọi là “overlap”. Việc xếp chồng các frame với nhau sẽ được thực hiện theo một tỷ lệ chồng lấp thích hợp, thông thường là 40% hoặc 50%. Sau khi các frame tín hiệu được xử lý triệt nhiễu trong miền tần số, các frame này được liên kết lại nhau bằng phương pháp thích hợp với phương pháp phân tích tín hiệu thành các frame ở đầu vào gọi là “adding”. Tập hợp các mẫu tín hiệu trong cùng một frame sau khi được phân tích ở đầu vào gọi là một segment. Với cách thực hiện phân tích và liên kết các frame bằng phương pháp overlap và adding thì tín hiệu của chúng ta thu được sau khi xử lý triệt nhiễu sẽ không bị méo dạng và sẽ không xuất hiện hiện tượng “giả nhiễu”. Hình 3.5 quá trình thực hiện overlap và adding [32]. Ước lượng và cập nhật nhiễu Phương thức ước lượng nhiễu có thể ảnh hưởng lớn đến chất lượng của tín hiệu sau khi được tăng cường. Nếu nhiễu được ước lượng quá nhỏ thì nhiễu sẽ vẫn còn trong tín hiệu và nó sẽ được nghe thấy, còn nếu như nhiễu được ước lượng quá lớn thì tiếng nói sẽ bị méo, và làm sẽ làm tính dễ nghe của tiếng nói bị ảnh hưởng. Cách đơn giản nhất để ước lượng và cập nhật phổ của nhiễu trong đoạn tín hiệu không có mặt của tiếng nói sử dụng thuật toán thăm dò hoạt động của tiếng nói (voice activity detection - VAD). Tuy nhiên phương pháp đó chỉ thoả mãn đối với nhiễu không thay đổi(nhiễu trắng), nó sẽ không hiệu quả trong các môi trường thực tế (ví dụ như nhà hàng), ở những nơi đó đặc tính phổ của nhiễu thay đổi liên tục. Trong mục này chúng ta sẽ đề cập đến thuật toán ước lượng nhiễu thay đổi liên tục và thực hiện trong lúc tiếng nói hoạt động, thuật toán này sẽ phù hợp môi trường có nhiễu thay đổi cao. Thuật toán thăm dò hoạt động của tiếng nói (VAD) Quá trình xử lý để phân biệt khi nào có tiếng nói hoạt động, khi nào không có tiếng nói (im lặng) được gọi là sự thăm dò hoạt động của tiếng nói – Voice activity detection (VAD). Thuật toán VAD có tín hiệu ra ở dạng nhị phân quyết định trên một nền tảng frame-by-frame, khi đó frame có thể xấp xỉ 20-40 ms. Một đoạn tiếng nói có chứa tiếng nói hoạt động thì VAD = 1, còn nếu tiếng nói không hoạt động hay đó chính là nhiễu thì VAD = 0. Có một vài thuật toán VAD được đưa ra dựa trên nhiều đặc tính của tín hiệu. Các thuật toán VAD được đưa ra sớm nhất thì dựa vào các đặc tính như mức năng lượng, zero-crossing, đặc tính cepstral, phép đo khoảng cách phổ Itakura LPC, phép đo chu kỳ. Phần lớn các thuật toán VAD đều phải đối mặt với vấn đề là điều kiện SNR thấp, đặc biệt khi nhiễu bị thay đổi. Một thuật toán VAD có độ chính xác trong môi trường thay đổi không thể đủ trong các ứng dụng của Speech enhancement, nhưng việc ước lượng nhiễu một cách chính xác là rất cần thiết tại mọi thời điểm khi tiếng nói hoạt động . Ước lượng và cập nhật nhiễu Nhiễu sẽ được ước lượng lúc ban đầu bằng cách lấy trung bình biên độ phổ của tín hiệu bị nhiễu : (3.28) Sau đó, sử dụng phương pháp VAD để nhận biết các frame tiếp theo, frame nào là frame nhiễu và sẽ cập nhật nhiễu đó cho các frame tiếp theo. Để có thể nhận biết được frame nào là nhiễu thì chúng ta thực hiện so sánh biên độ phổ của nhiễu được ước lượng với biên độ phổ của tín hiệu bị nhiễu : (3.29) Nếu thì frame đó không phải là frame có tiếng nói, khi đó ta có thể cập nhật lại nhiễu đã được ước lượng trước đó. Kết luận Để thuật toán Wiener Filter có thể thực hiện được thì cần phải phân tích tín hiệu thành các frame và các frame phải xếp chồng lên nhau, và sau khi các frame được xử lý trong miền tần số và chuyển đổi về lại miền thời gian thì các frame đó phải được liên kết lại với nhau theo đúng phương pháp tương ứng với phương pháp phân tích tín hiệu ở đầu vào, quá trình đó gọi là overlap và adding. Chính điều đó sẽ làm cho tín hiệu của chúng ta sau khi xử lý triệt nhiễu sẽ không bị méo, đảm bảo chất lượng của tiếng nói. Nội dung của chương cũng trình bày vấn đề ước lượng nhiễu, đây là cái chính mà speech enhancement cần giải quyết, nó quyết định tính hiệu quả của thuật toán và chất lượng của tiếng nói sau khi xử lý triệt nhiễu. Thực hiện thuật toán Wiener trên Matlab Dựa vào lý thuyết đã nghiên cứu phần trên,phần này chúng ta xây dựng các sơ đồ thuật toán và thực hiện các thuật toán giảm nhiễu bằng Matlab. Sơ đồ thuật toán Wiener Filter Tính lại mức nhiễu N End I=I+1;nhập frame tiếp theo Begin Phân chia Frame tín hiệu đầu vào Tinh cong suat nhieu trung binh N ban đầu I=0;Nhập frame đầu tiên VAD X(:,i)=Beta*Y(:,i) D=YS(:,i)-N; % Thực hiện trừ phổ X(:,i)=max(D,0); Y=biến đổi FFT cho các frame X = X = Đ = X = S Đ SpeechFlag==0? S I<number of frame Thực hiên IFFT và nối các frame Đ Hình 2 : Sơ đồ thuật toán Wiener Filter Thực hiện bằng Matlab Chương trình function output=WienerScalart96(signal,fs,IS) % output=WIENERSCALART96(signal,fs,IS) % Wiener filter based on tracking a priori SNR usingDecision-Directed % method, proposed by Scalart et al 96. In this method it is assumed that % SNRpost=SNRprior +1. based on this the Wiener Filter can be adapted to a % model like Ephraims model in which we have a gain function which is a % function of a priori SNR and a priori SNR is being tracked using Decision % Directed method. % Author: Esfandiar Zavarehei % Created: MAR-05 if (nargin<3 | isstruct(IS)) IS=.25; %Initial Silence or Noise Only part in seconds end W=fix(.025*fs); %Window length is 25 ms SP=.4; %Shift percentage is 40% (10ms) %Overlap-Add method works good with this value(.4) wnd=hamming(W); %IGNORE FROM HERE ............................... if (nargin>=3 & isstruct(IS))%This option is for compatibility with another programme W=IS.windowsize SP=IS.shiftsize/W; %nfft=IS.nfft; wnd=IS.window; if isfield(IS,'IS') IS=IS.IS; else IS=.25; end end % ......................................UP TO HERE pre_emph=0; signal=filter([1 -pre_emph],1,signal); NIS=fix((IS*fs-W)/(SP*W) +1);%number of initial silence segments y=segment(signal,W,SP,wnd); % This function chops the signal into frames Y=fft(y); YPhase=angle(Y(1:fix(end/2)+1,:)); %Noisy Speech Phase Y=abs(Y(1:fix(end/2)+1,:));%Specrogram numberOfFrames=size(Y,2); FreqResol=size(Y,1); N=mean(Y(:,1:NIS)')'; %initial Noise Power Spectrum mean LambdaD=mean((Y(:,1:NIS)').^2)';%initial Noise Power Spectrum variance alpha=.99; %used in smoothing xi (For Deciesion Directed method for estimation of A Priori SNR) NoiseCounter=0; NoiseLength=9;%This is a smoothing factor for the noise updating G=ones(size(N));%Initial Gain used in calculation of the new xi Gamma=G; X=zeros(size(Y)); % Initialize X (memory allocation) h=waitbar(0,'Wait...'); for i=1:numberOfFrames %%%%%%%%%%%%%%%%VAD and Noise Estimation START if i<=NIS % If initial silence ignore VAD SpeechFlag=0; NoiseCounter=100; else % Else Do VAD [NoiseFlag, SpeechFlag, NoiseCounter, Dist]=vad(Y(:,i),N,NoiseCounter); %Magnitude Spectrum Distance VAD end if SpeechFlag==0 % If not Speech Update Noise Parameters N=(NoiseLength*N+Y(:,i))/(NoiseLength+1); %Update and smooth noise mean LambdaD=(NoiseLength*LambdaD+(Y(:,i).^2))./(1+NoiseLength); %Update and smooth noise variance end %%%%%%%%%%%%%%%%%%%VAD and Noise Estimation END gammaNew=(Y(:,i).^2)./LambdaD; %A postiriori SNR xi=alpha*(G.^2).*Gamma+(1-alpha).*max(gammaNew-1,0); %Decision Directed Method for A Priori SNR Gamma=gammaNew; G=(xi./(xi+1)); X(:,i)=G.*Y(:,i); %Obtain the new Cleaned value waitbar(i/numberOfFrames,h,num2str(fix(100*i/numberOfFrames))); end close(h); output=OverlapAdd2(X,YPhase,W,SP*W); %Overlap-add Synthesis of speech output=filter(1,[1 -pre_emph],output); %Undo the effect of Pre-emphasis output=0.999*(output/max(abs(output))); function ReconstructedSignal=OverlapAdd2(XNEW,yphase,windowLen,ShiftLen); %Y=OverlapAdd(X,A,W,S); %Y is the signal reconstructed signal from its spectrogram. X is a matrix %with each column being the fft of a segment of signal. A is the phase %angle of the spectrum which should have the same dimension as X. if it is %not given the phase angle of X is used which in the case of real values is %zero (assuming that its the magnitude). W is the window length of time %domain segments if not given the length is assumed to be twice as long as %fft window length. S is the shift length of the segmentation process ( for %example in the case of non overlapping signals it is equal to W and in the %case of %50 overlap is equal to W/2. if not givven W/2 is used. Y is the %reconstructed time domain signal. %Sep-04 %Esfandiar Zavarehei if nargin<2 yphase=angle(XNEW); end if nargin<3 windowLen=size(XNEW,1)*2; end if nargin<4 ShiftLen=windowLen/2; end if fix(ShiftLen)~=ShiftLen ShiftLen=fix(ShiftLen); disp('The shift length have to be an integer as it is the number of samples.') disp(['shift length is fixed to ' num2str(ShiftLen)]) end [FreqRes FrameNum]=size(XNEW); Spec=XNEW.*exp(j*yphase); if mod(windowLen,2) %if FreqResol is odd Spec=[Spec;flipud(conj(Spec(2:end,:)))]; else Spec=[Spec;flipud(conj(Spec(2:end-1,:)))]; end sig=zeros((FrameNum-1)*ShiftLen+windowLen,1); weight=sig; for i=1:FrameNum start=(i-1)*ShiftLen+1; spec=Spec(:,i); sig(start:start+windowLen-1)=sig(start:start+windowLen-1)+real(ifft(spec,windowLen)); end ReconstructedSignal=sig; function Seg=segment(signal,W,SP,Window) % SEGMENT chops a signal to overlapping windowed segments % A= SEGMENT(X,W,SP,WIN) returns a matrix which its columns are segmented % and windowed frames of the input one dimentional signal, X. W is the % number of samples per window, default value W=256. SP is the shift % percentage, default value SP=0.4. WIN is the window that is multiplied by % each segment and its length should be W. the default window is hamming % window. % 06-Sep-04 % Esfandiar Zavarehei if nargin<3 SP=.4; end if nargin<2 W=256; end if nargin<4 Window=hamming(W); end Window=Window(:); %make it a column vector L=length(signal); SP=fix(W.*SP); N=fix((L-W)/SP +1); %number of segments Index=(repmat(1:W,N,1)+repmat((0:(N-1))'*SP,1,W))'; hw=repmat(Window,1,N); Seg=signal(Index).*hw; function [NoiseFlag, SpeechFlag, NoiseCounter, Dist]=vad(signal,noise,NoiseCounter,NoiseMargin,Hangover) %[NOISEFLAG, SPEECHFLAG, NOISECOUNTER, DIST]=vad(SIGNAL,NOISE,NOISECOUNTER,NOISEMARGIN,HANGOVER) %Spectral Distance Voice Activity Detector %SIGNAL is the the current frames magnitude spectrum which is to labeld as %noise or speech, NOISE is noise magnitude spectrum template (estimation), %NOISECOUNTER is the number of imediate previous noise frames, NOISEMARGIN %(default 3)is the spectral distance threshold. HANGOVER ( default 8 )is %the number of noise segments after which the SPEECHFLAG is reset (goes to %zero). NOISEFLAG is set to one if the the segment is labeld as noise %NOISECOUNTER returns the number of previous noise segments, this value is %reset (to zero) whenever a speech segment is detected. DIST is the %spectral distance. %Saeed Vaseghi %edited by Esfandiar Zavarehei %Sep-04 if nargin<4 NoiseMargin=3; end if nargin<5 Hangover=8; end if nargin<3 NoiseCounter=0; end FreqResol=length(signal); SpectralDist= 20*(log10(signal)-log10(noise)); SpectralDist(find(SpectralDist<0))=0; Dist=mean(SpectralDist); if (Dist < NoiseMargin) NoiseFlag=1; NoiseCounter=NoiseCounter+1; else NoiseFlag=0; NoiseCounter=0; end % Detect noise only periods and attenuate the signal if (NoiseCounter > Hangover) SpeechFlag=0; else SpeechFlag=1; end Giải thích chương trình theo thuật toán function output=WienerScalart96(signal,fs,IS) Tạo hàm WienerScalart96 với 3 tham số đầu vào là tín hiệu signal ,tần số lấy mẫu fs và một tham số IS. Phân chia Frame tín hiệu đầu