extract a thin band of frequencies from entire spectrum using FFT

Hello,

I have following program:

fs=1e4;
x=randn(1,1e4); %gaussian dist random signal
x_fft=fft(x); %fft of the signal
N=length(x);
omega=fs*[(0:N/2) (-N/2+1:-1)]/N; %frequency axis(bins)
figure;
plot(omega(1:N/2),2*abs(x_fft(1:N/2))); %plot of one-sided spectrum

Now I want to extract only a thin band of frequencies(say 500-520Hz) from the spectrum,take its ifft store the obtained time-domain signal in a variable y. I know it is bandpass filtering in TD and there is function fdesign.bandpass() in matlab aswell for this purpose, but I want to do it in frequency domain (using FFT and IFFT) as in my case, the samples are very large and BPF is very time consuming.

At the end, I want to have a small function, where I only enter the frequency band of interest as my input and I get the corresponding TD signal extracted as my output.

Please help.

Answer 1

[Azzi Abdelmalek]

fs=1e4;pas=1/fs
x=randn(1,1e4); %gaussian dist random signal
x_fft=fft(x); %fft of the signal
N=length(x);
[wc,w0,a0,ak,bk,c0,ck]=get_harmonics(x_fft,pas,2)
fc=wc/(2*pi),

%to extract index frequencies from ck

idx=find(fc>=500 & fc<=520)

% corresponding temporel signal signal(t)

t=0:0.1:10000;
signal=zeros(1,numel(t));
for l=1:length(idx)
  k=idx(l);
  signal=signal+ck(k)*exp(j*k*w0*t);
end
plot(t,real(signal)); %more t is near inf more imag(signal) is near 0.

% download get_harmonics at http://www.mathworks.com/matlabcentral/fileexchange/37654-get-harmoniques-of-a-real-signal

[zozo]

Thank you @Azzi, but When I plot:

plot(abs(fft(signal)))

it is not conjugate symmetric..why? it is supposed to be conjugate right?

[Azzi Abdelmalek]

The fourier tranform of a real continuous signal is conjugate symetric, the FFT is used for discret signals, and when we use it for continuous signals, we have to make some modifications to obtain our result

Answer 2

[Wayne King]

You can just create a vector of zeros to match the size of the original DFT vector and then fill the correct indices of the vectors with the DFT coefficients. Remember for a real-valued signal, you will have to fill the appropriate negative and positive frequencies with the complex conjugates. Here I'll just give a simple example to pull out the 100-Hz component. For a band of frequencies, you'll need the vector for the band.

 t = 0:0.001:1-0.001;
 x = cos(2*pi*100*t)+cos(2*pi*200*t);
 xdft = fft(x);
 % the DFT bins for 100 Hz are 101 and 1000-101+2
 indices100 = [101 length(xdft)-101+2];
 ydft = zeros(size(xdft));
 ydft(indices100) = xdft(indices100);
 y = ifft(ydft,'symmetric');

[zozo]

Thank you @Wayne