Alguma regra prudente como suavizar o espectro da FFT para evitar artefatos ao fazer a mão?

StackOverflow https://stackoverflow.com/questions/3827346

  •  26-09-2019
  •  | 
  •  

Pergunta

Eu tenho um espectro de magnitude da FFT e quero criar um filtro a partir dele que passa seletivamente fontes de ruído periódicas (por exemplo, esporas de ondas tendrâneas) e o Zero está fora das caixas de frequência associadas ao ruído de fundo aleatório. Entendo que transições nítidas no domínio FREQ criarão artefatos de toque quando esse filtro estiver de volta ao domínio do tempo ... e por isso estou me perguntando se há alguma regras práticas como suavizar as transições em tal filtro para evitar tais tais toque.

Por exemplo, se a FFT tiver 1M de caixas de frequência e há cinco esporões saindo do piso de ruído de fundo, gostaria de zerar todas as caixas, exceto a lixeira de pico associada a cada um dos cinco esporões. A questão é como lidar com as caixas de esporão vizinhas para evitar artefatos no domínio do tempo. Por exemplo, a lixeira de cada lado de uma lixeira deve ser ajustada para 50% de amplitude? Devem ser usadas duas caixas em ambos os lados de uma lixeira (a mais próxima a 50%e a próxima mais próxima a 25%, etc.)? Quaisquer pensamentos muito apreciados. Obrigado!

Foi útil?

Solução

Eu gosto do seguinte método:

  • Crie o espectro de magnitude ideal (lembrando -se de torná -lo simétrico sobre DC)
  • Transformação inversa no domínio do tempo
  • Gire o bloco pela metade do tamanho do bloco
  • Aplique uma janela Hann

Acho que cria resultados de domínio de frequência razoavelmente suave, embora eu nunca tenha experimentado algo tão nítido quanto você está sugerindo. Você provavelmente pode fazer um filtro mais nítido usando uma janela Kaiser-Bessel, mas precisa escolher os parâmetros adequadamente. Por Sharper, acho que talvez você possa reduzir os lóbulos laterais em 6 dB ou mais.

Aqui estão alguns amostras de código MATLAB/OCTAVE. Para testar os resultados, usei freqz(h, 1, length(h)*10);.

function [ht, htrot, htwin] = ArbBandPass(N, freqs)
%# N = desired filter length
%# freqs = array of frequencies, normalized by pi, to turn into passbands
%# returns raw, rotated, and rotated+windowed coeffs in time domain

if any(freqs >= 1) || any(freqs <= 0)
    error('0 < passband frequency < 1.0 required to fit within (DC,pi)')
end

hf = zeros(N,1); %# magnitude spectrum from DC to 2*pi is intialized to 0
%# In Matlabs FFT, idx 1 -> DC, idx 2 -> bin 1, idx N/2 -> Fs/2 - 1, idx N/2 + 1 -> Fs/2, idx N -> bin -1
idxs = round(freqs * N/2)+1; %# indeces of passband freqs between DC and pi
hf(idxs) = 1; %# set desired positive frequencies to 1
hf(N - (idxs-2)) = 1; %# make sure 2-sided spectrum is symmetric, guarantees real filter coeffs in time domain
ht = ifft(hf); %# this will have a small imaginary part due to numerical error
if any(abs(imag(ht)) > 2*eps(max(abs(real(ht)))))
    warning('Imaginary part of time domain signal surprisingly large - is the spectrum symmetric?')
end
ht = real(ht); %# discard tiny imag part from numerical error
htrot = [ht((N/2 + 1):end) ; ht(1:(N/2))]; %# circularly rotate time domain block by N/2 points
win = hann(N, 'periodic'); %# might want to use a window with a flatter mainlobe
htwin = htrot .* win;
htwin = htwin .* (N/sum(win)); %# normalize peak amplitude by compensating for width of window lineshape
Licenciado em: CC-BY-SA com atribuição
Não afiliado a StackOverflow
scroll top