www.pudn.com > MATLAB.rar > noise_inv_f.m, change:2010-01-07,size:2066b

function hout=noise_inv_f(L); % L= filter length/2 % Create a filter kernel (FIR coefficients) that approximates a 1/f frequency response. % Dick Benson September 2004 % Copywrite 2004-2010 The MathWorks, Inc. if nargin==0 L = 128; % filter kernel length/2 end; f =(0:(L-1))/L; % frequency vector normalized to Fs/2 shape = (1./(f+0.1/L)).^0.5; % desired shape of spectrum (note +0.1/L to limit infinity ) , % added .^0.5 to make POWER drop as 1/f shape = shape/shape(1); % normalize DC gain to 1 (0dB) ph = pi*(0:L-1); % phase shift to center impulse response in time window % Need to create a complex spectral description % which, when transformed with the ifft, % creates a REAL impulse response. mag=[shape,0, fliplr(shape(2:end))]; % construct magnitude phase=[-ph,0,fliplr(ph(2:end))]; % and phase H=mag.*exp(j*phase); % complex representation h=ifft(H); % impulse reponse % sanity checks .... subplot(2,1,1) plot((h)); title('1/f filter impulse response (no window)'); xlabel('samples') % If the filter is used as is, there will be ripples in the frequency reponse, % due to the abrupt transistions at the impulse response endpoints. % To mitigate this, apply a Hanning window and then compare final response with % desired. hout=h.*hann(2*L)'; H2=fft(hout); % take fft of windowed filter kernel % hf=figure; subplot(2,1,2) semilogx(f/2,20*log10(abs(H2(1:L))),f/2,20*log10(shape),'x'); legend('attained','desired') title('overlayed freq responses of desired and attained 1/f shape'); xlabel('Normalized to Fs'); ylabel('dB') % Finally, scale the filter so that 1 Watt of (broad band) input noise % produces 1 Watt of output noise. df = f(2)-f(1); Pout_T = sum(hout.^2); % power out for 1 W input hout = hout/sqrt(Pout_T); % use filter_check.mdl to confirm power level and shape f_inv = hout; save f_inv f_inv; % save for use by other examples. set(gcf,'HandleVisibility','off');