Thread Subject:

FFT over a limited frequency range

Hi,

I want to analyze the first 5 Hz of a signal spectrum..I've used the fft to calculate the spectrum, but it is 70 Hz wide, so I have a lot of "useless" frequency..

It's possible to "concentrate" the 1024 FFT points (number I'm actually using) in the first (about) 5 Hz?
In this way I could get an higher frequency resolution, higher than the actual 0.1 Hz.

I've tried the Goertzel algorithm, but it works only with integer frequency indexes, so it doesn't work fine for me..

Any idea?

Thanks

"Mario Rossi" <nonono@hotmail.com> wrote in message <idqg3a$ljn$1@ginger.mathworks.com>...
> Hi,
>
> I want to analyze the first 5 Hz of a signal spectrum..I've used the fft to calculate the spectrum, but it is 70 Hz wide, so I have a lot of "useless" frequency..
>
> It's possible to "concentrate" the 1024 FFT points (number I'm actually using) in the first (about) 5 Hz?
> In this way I could get an higher frequency resolution, higher than the actual 0.1 Hz.
>
> I've tried the Goertzel algorithm, but it works only with integer frequency indexes, so it doesn't work fine for me..
>
> Any idea?

Hi Mario,

You could try to filter out any frequency components above the first 5Hz band. FFT will always give you a decomposition of whatever is present in the signal, so if those higher frequency components are there, they wont show up in the FFT.

Hope this helps.

Nic
>
> Thanks

"Mario Rossi" <nonono@hotmail.com> wrote in message <idqg3a$ljn$1@ginger.mathworks.com>...
> Hi,
>
> I want to analyze the first 5 Hz of a signal spectrum..I've used the fft to calculate the spectrum, but it is 70 Hz wide, so I have a lot of "useless" frequency..
>
> It's possible to "concentrate" the 1024 FFT points (number I'm actually using) in the first (about) 5 Hz?
> In this way I could get an higher frequency resolution, higher than the actual 0.1 Hz.
>
> I've tried the Goertzel algorithm, but it works only with integer frequency indexes, so it doesn't work fine for me..
>
> Any idea?
>
> Thanks

Hi Mario, the frequency bins (in Hz) in the DFT are determined by the number of samples and the sampling frequency, so you can't improve your resolution in the way you are attempting.

However, if you have the Signal Processing Toolbox, there are superresolution methods (subspace methods) for frequency estimation, which may serve your purpose.

See spectrum.music for example.

Wayne

"Nic Roberts" <dingtheking@googlemail.com> wrote in message <idqkfg$2j2$1@fred.mathworks.com>...
>
> You could try to filter out any frequency components above the first 5Hz band. FFT will always give you a decomposition of whatever is present in the signal, so if those higher frequency components are there, they wont show up in the FFT.
======

And then, I'm guessing, you would redo the FFT with a larger time sampling interval, dT, in order to cover only the smaller sampled bandwidth of 5 Hz. That's because of the relationship

dT=1/sampledBandwidth

Hi guys,

thank you for your advices..

I'm searching for some documentation on superresolution methods..then I'll post something if I found a solution..

For the other suggestion: I can't change my sampling frequency, since it's hardware depending... If I use zero-padding it could be useful or not?

Thanks

"Mario Rossi" <nonono@hotmail.com> wrote in message <ie30l9$ir$1@fred.mathworks.com>...
> Hi guys,
>
> thank you for your advices..
>
> I'm searching for some documentation on superresolution methods..then I'll post something if I found a solution..
>
> For the other suggestion: I can't change my sampling frequency, since it's hardware depending... If I use zero-padding it could be useful or not?
>
> Thanks

Hi Mario,

Zero padding (in the time domain) would show up as interpolation in the freq domain, so in a sense it could "increase" the resolution between freq samples, but is in reality just interpolation, it could be useful sometimes.

Regards,
Juan

Mario Rossi wrote:
> Hi guys,
>
> thank you for your advices..
>
> I'm searching for some documentation on superresolution methods..then
> I'll post something if I found a solution..
>
> For the other suggestion: I can't change my sampling frequency, since
> it's hardware depending... If I use zero-padding it could be useful or not?
...

Also look at "Zoom" FFT...here's a quick overview link (aimed at speech
processing, but concept applies regardless of application/bandwidth)

<http://www.numerix-dsp.com/zoomfft.html>

W/ baseband FFT and a constrained sample rate, other than the
interpolation technique of zero-padding, you _could_ increase the sample
time as the resulting df = 1/T. This is more real information at the
cost of a longer time series and, of course, longer FFT as well (covers
same overall Nyquist, just finer resolution).

--

Hi,

thank you again, I'm trying with your suggestion but I have no solutions yet.
So far, I've tried:

- Goertzel algorithm: it seems perfect, but it works only on integer frequency... isn't it? At least the matlab implementation accept only integer values, and I need about 0.01 Hz of resolution..

- ZoomFFT: I've found a lot of implementation, but it doesn't seems to be so useful.

- I'm trying this method:
calculate the dft directly, without any algorithm:

f1=0.01;
f2=10;
step=0.01;
spectrum = zeros (1,uint32(((f2-f1)/step)+1));
frequencyRange=(f1:passo:f2);

% this is the signal to analyze
vector = filtrato (:,343)';

% classical DFT calculation
for f=f1:step:f2
       for n=1:length(vector)
            spectrum(1,f/step) = spectrum(1,f/step) + vector(n)*exp(-1i*2*pi*n*f*step));
       end;
end;

plot(frequencyRange,abs(spectrum));

It isn't so bad... but maybe I make some mistake... here are the plots, with this manual DFT and the classic matlab FFT:
[URL=http://img27.imageshack.us/i/dfty.jpg/][IMG]http://img27.imageshack.us/img27/3931/dfty.jpg[/IMG][/URL]
[URL=http://img405.imageshack.us/i/fftx.jpg/][IMG]http://img405.imageshack.us/img405/2315/fftx.jpg[/IMG][/URL]

the frequency are wrong...
And even worse if I increase the resolution (and so if I reduce the step):
[URL=http://img412.imageshack.us/i/dftmoreresolution.jpg/][IMG]http://img412.imageshack.us/img412/7030/dftmoreresolution.jpg[/IMG][/URL]


Any idea?
thank you

Sorry, links don't work.. Here they are:

DFT and FFT
[URL=http://img27.imageshack.us/i/dfty.jpg/[/URL]
[URL=http://img405.imageshack.us/i/fftx.jpg/[/URL]

More resolution DFT
[URL=http://img412.imageshack.us/i/dftmoreresolution.jpg/[/URL]

They should work now..

On Dec 30, 5:30 am, "Mario Rossi" <non...@hotmail.com> wrote:
> Hi,
>
> thank you again, I'm trying with your suggestion but I have no solutions yet.
> So far, I've tried:
>
> - Goertzel algorithm: it seems perfect, but it works only on integer frequency... isn't it? At least the matlab implementation accept only integer values, and I need about 0.01 Hz of resolution..
>
> - ZoomFFT: I've found a lot of implementation, but it doesn't seems to be so useful.
>
> - I'm trying this method:
> calculate the dft directly, without any algorithm:
>
> f1=0.01;
> f2=10;
> step=0.01;
> spectrum = zeros (1,uint32(((f2-f1)/step)+1));
> frequencyRange=(f1:passo:f2);
>
> % this is the signal to analyze
> vector = filtrato (:,343)';
>
> % classical DFT calculation
> for f=f1:step:f2
>        for n=1:length(vector)
>             spectrum(1,f/step) = spectrum(1,f/step) + vector(n)*exp(-1i*2*pi*n*f*step));    
>        end;
> end;
>
> plot(frequencyRange,abs(spectrum));
>
> It isn't so bad... but maybe I make some mistake... here are the plots, with this manual DFT and the classic matlab FFT:
> [URL=http://img27.imageshack.us/i/dfty.jpg/][IMG]http://img27.imageshack.us/img27/3931/dfty.jpg[/IMG][/URL]
> [URL=http://img405.imageshack.us/i/fftx.jpg/][IMG]http://img405.imageshack.us/img405/2315/fftx.jpg[/IMG][/URL]
>
> the frequency are wrong...
> And even worse if I increase the resolution (and so if I reduce the step):
> [URL=http://img412.imageshack.us/i/dftmoreresolution.jpg/][IMG]http://img412.imageshack.us/img412/7030/dftmoreresolution.jpg[/IMG][/URL]
>
> Any idea?
> thank you

http://groups.google.com/group/comp.soft-sys.matlab/msg/be4844a338a5f500

function [Xdft,Xqr,Xpi,NMSEx] = dftgh6(x,t,f)

% function [Xdft,Xqr,Xpi,NMSEx] = dftgh6(x,t,f)
%
% Extension of AJ Johnson's dft to NONUNIFORM sampling
%
% Compute Xdft (Discrete Fourier Transform) at M frequencies
% given in f, given N samples x taken at times t:
%
% Xdft(f) = sum(k=1,N){dts(k)*x(k)*exp(-j*2*pi*f*t(k)*f)}
% = W *(x.*dts)
% where dts is a symmetrized modification of diff(t).
%
% Compute least squares spectra Xqr and Xpi using
% using QR backslash and pseudoinverse least square
% inversion, respectively.
%
% NMSEx is the normalized mean-square-error of
% reconstucting x from the Xdft, Xqr and Xpi spectra.
% If meanx = x, i.e., x is constant, then the MSE is
% unnormalized.
%
% NMSEx(1) error from Xqr => xidftqr
% NMSEx(2) error from Xpi => xidftpi
% NMSEx(3) error from Xdft => xiqrdft
% NMSEx(4) error from Xdft => xipidft
%
% For comparison with MATLAB's Xfft when the spacing is
% uniform, multiply Xfft by dt0 = mean(diff(t))
%
% Greg Heath 15 July 2010

x = x(:); % Format 'x' into a column vector
t = t(:); % Format 't' into a column vector
f = f(:); % Format 'f' into a column vector
N = length(x);
if length(t)~= N
   error('x and t do not have the same length')
end;

dt = diff(t); % asymmetric "dt"
% symmetric "dt"
dts = [dt(1);(dt(1:end-1)+dt(2:end))/2;dt(end)];
meanx = sum(x.*dts)/sum(dts);

W = exp(-2*pi*j * f*t');
Xdft = W * (x.*dts);

df = diff(f); % asymmetric "df"
% symmetric "df"
dfs = [df(1);(df(1:end-1)+df(2:end))/2;df(end)];

% x = (1/N) *W' *(X.*dfs) IDFT

Xqr = ( W'\x )./dfs;
Xpi = (pinv(W')*x)./dfs;

xidftqr = W' *(Xqr.*dfs); % Fourier Reconstruction
xidftpi = W' *(Xpi.*dfs); % Fourier Reconstruction

xiqrdft = ( W\Xdft )./dts ; % QR LSE Reconstruction
% xipidft = ( pinv(W)*Xdft )./dts ; % PI LSE Reconstruction

MSE0 = mse(abs(x-meanx));
if MSE0 == 0, MSE0 = 1, end;

NMSEx(1)= mse(abs(x-xidftqr))/MSE0;
NMSEx(2)= mse(abs(x-xidftpi))/MSE0;
NMSEx(3)= mse(abs(x-xiqrdft))/MSE0;
NMSEx(4)= mse(abs(x-xipidft))/MSE0;

return

Hope this helps.

Greg

On Dec 9, 6:55 am, "Mario Rossi" <non...@hotmail.com> wrote:
> Hi,
>
> I want to analyze the first 5 Hz of a signal spectrum..I've used the fft to calculate the spectrum, but it is 70 Hz wide, so I have a lot of "useless" frequency..
>
> It's possible to "concentrate" the 1024 FFT points (number I'm actually using) in the first (about) 5 Hz?
> In this way I could get an higher frequency resolution, higher than the actual 0.1 Hz.

I'm confused.

To get a resolution of df = 0.1 Hz with N = 1024
points, you need a sampling frequency of

Fs = N*df = 102.4
dt = 1/Fs = 9.766 msec
tmax = (N-1)*dt = 9.99 sec
fmax = Fs/2 = 51.2 Hz

Where does "70 Hz" come from?

> I've tried the Goertzel algorithm, but it works
> only with integer frequency indexes, so it doesn't
> work fine for me..

Almost all DSP problems involve integer frequency indexing.

Why is this a problem?

Greg

Greg Heath <heath@alumni.brown.edu> wrote in message <94453d96-2561-4fdf-aad7-9bb451170366@c2g2000yqc.googlegroups.com>...
> On Dec 9, 6:55 am, "Mario Rossi" <non...@hotmail.com> wrote:
> > Hi,
> >
> > I want to analyze the first 5 Hz of a signal spectrum..I've used the fft to calculate the spectrum, but it is 70 Hz wide, so I have a lot of "useless" frequency..
> >
> > It's possible to "concentrate" the 1024 FFT points (number I'm actually using) in the first (about) 5 Hz?
> > In this way I could get an higher frequency resolution, higher than the actual 0.1 Hz.
>
> I'm confused.
>
> To get a resolution of df = 0.1 Hz with N = 1024
> points, you need a sampling frequency of
>
> Fs = N*df = 102.4
> dt = 1/Fs = 9.766 msec
> tmax = (N-1)*dt = 9.99 sec
> fmax = Fs/2 = 51.2 Hz
>
> Where does "70 Hz" come from?
>
> > I've tried the Goertzel algorithm, but it works
> > only with integer frequency indexes, so it doesn't
> > work fine for me..
>
> Almost all DSP problems involve integer frequency indexing.
>
> Why is this a problem?
>
> Greg

Hi Greg,
thank you for your code, that I'm going to analyze in these days..

About your second post..
70 Hz come from the length of the signal = 14 s, and the number of samples (1000) in this signal.
fs = 1/(14/1000) = 71 Hz (more or less)
I was wrong when I said "actual resolution of 0.1 Hz", it is about 0.07 Hz.
However the problem is still the same, I have to improve the resolution..

Goertzel:
for example, I'd like to use Goertzel, in a for cicle, to calculate all the single spectrum components in a range from 0,01 Hz and 5 Hz.
But if I pass to the algorithm the non integer values (all values apart from 1, 2, 3, 4 ,5 Hz!!!) it returns an error! It says that it accepts only integer values..
Do you know a way to avoid this problem?

Thanks,

Mario

On Dec 31, 4:25 am, "Mario Rossi" <non...@hotmail.com> wrote:
> Greg Heath <he...@alumni.brown.edu> wrote in message <94453d96-2561-4fdf-aad7-9bb451170...@c2g2000yqc.googlegroups.com>...
> > On Dec 9, 6:55 am, "Mario Rossi" <non...@hotmail.com> wrote:
> > > Hi,
>
> > > I want to analyze the first 5 Hz of a signal spectrum..I've used the fft to calculate the spectrum, but it is 70 Hz wide, so I have a lot of "useless" frequency..
>
> > > It's possible to "concentrate" the 1024 FFT points (number I'm actually using) in the first (about) 5 Hz?
> > > In this way I could get an higher frequency resolution, higher than the actual 0.1 Hz.
>
> > I'm confused.
>
> > To get a resolution of df = 0.1 Hz with N = 1024
> > points, you need a sampling frequency of
>
> > Fs   = N*df     = 102.4
> > dt   = 1/Fs     = 9.766 msec
> > tmax = (N-1)*dt = 9.99 sec
> > fmax = Fs/2     = 51.2 Hz
>
> > Where does "70 Hz" come from?
>
> > > I've tried the Goertzel algorithm, but it works
> > > only with integer frequency indexes, so it doesn't
> > > work fine for me..
>
> > Almost all DSP problems involve integer frequency indexing.
>
> > Why is this a problem?
>
> > Greg
>
> Hi Greg,
> thank you for your code, that I'm going to analyze in these days..
>
> About your second post..
> 70 Hz come from the length of the signal = 14 s, and the number of samples (1000) in this signal.
> fs = 1/(14/1000) = 71 Hz (more or less)
> I was wrong when I said "actual resolution of 0.1 Hz", it is about 0.07 Hz.
> However the problem is still the same, I have to improve the resolution..
>
> Goertzel:
> for example, I'd like to use Goertzel, in a for cicle, to calculate all the single spectrum components in a range from 0,01 Hz and 5 Hz.
> But if I pass to the algorithm the non integer values (all values apart from 1, 2, 3, 4 ,5 Hz!!!) it returns an error! It says that it accepts only integer values..
> Do you know a way to avoid this problem?
>
> Thanks,
>
> Mario- Hide quoted text -
>
> - Show quoted text -

Why not just use a dft? If df is constant you can use horners
method to obtain a much faster algorithm.

I will return tonight. Have to go look at a house.

Greg

On Dec 31 2010, 4:25 am, "Mario Rossi" <non...@hotmail.com> wrote:
-----SNIP
> Hi Greg,
> thank you for your code, that I'm going to analyze in these days..
>
> About your second post..
> 70 Hz come from the length of the signal = 14 s, and the number
> of samples (1000) in this signal.
> fs = 1/(14/1000) = 71 Hz (more or less)
> I was wrong when I said "actual resolution of 0.1 Hz", it is
> about 0.07 Hz.

Why approximate?

t = dt*(0:N-1);
tmax = dt*(N-1)
Fs = 1/dt = (N-1)/tmax = 71.357
df = 1/tmax = 0.071429

> However the problem is still the same, I have to improve the
> resolution..

Unless you use superresolution techniques like maximum entropy,
the only way to to that is to increase tmax.

> Goertzel:
> for example, I'd like to use Goertzel, in a for cicle, to
> calculate all the single spectrum components in a range from
> 0,01 Hz and 5 Hz.
> But if I pass to the algorithm the non integer values (all
> values apart from 1, 2, 3, 4 ,5 Hz!!!) it returns an error!
> It says that it accepts only integer values..
> Do you know a way to avoid this problem?

I just modified Maurice Givens' codes in

http://groups.google.com/group/comp.dsp/msg/14c9c255d1005396

Will post them after dinner.

Hope his helps.

Greg

Try the "Chirp Z Transform":

    doc czt


HTH.


Rick


"Mario Rossi" <nonono@hotmail.com> wrote in message
news:idqg3a$ljn$1@ginger.mathworks.com...
> Hi,
>
> I want to analyze the first 5 Hz of a signal spectrum..I've used the fft
> to calculate the spectrum, but it is 70 Hz wide, so I have a lot of
> "useless" frequency..
>
> It's possible to "concentrate" the 1024 FFT points (number I'm actually
> using) in the first (about) 5 Hz? In this way I could get an higher
> frequency resolution, higher than the actual 0.1 Hz.
>
> I've tried the Goertzel algorithm, but it works only with integer
> frequency indexes, so it doesn't work fine for me..
>
> Any idea?
> Thanks

On Jan 1, 6:20 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Dec 31 2010, 4:25 am, "Mario Rossi" <non...@hotmail.com> wrote:
> -----SNIP
> > Goertzel:
> > for example, I'd like to use Goertzel, in a for cicle, to
> > calculate all the single spectrum components in a range from
> > 0,01 Hz and 5 Hz.
> > But if I pass to the algorithm the non integer values (all
> > values apart from 1, 2, 3, 4 ,5 Hz!!!) it returns an error!
> > It says that it accepts only integer values..
> > Do you know a way to avoid this problem?
>
> I just modified Maurice Givens' codes in
>
> http://groups.google.com/group/comp.dsp/msg/14c9c255d1005396
>
> Will post them after dinner.

                  GOERTZELGH1

goertzelgh1 is used to detect a single frequency:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = goertzelgh1(x,Fs,f0)

% Modified code of Maurice Givens
% http://groups.google.com/group/comp.dsp/msg/...
% 14c9c255d1005396
%
% x = input vector
% Fs = input sampling frequency
% f0 = single test frequency
% A = test frequency amplplitude
%
% Greg Heath 1 Jan 2010

if f0 >= Fs/2
    error( 'f0 >= Nyquist frequency')
end

N = length(x);
coeff = 2 * cos(2*pi*f0/Fs);
Q1 = 0; Q2 = 0;
for j = 1:N
    Q0 = coeff * Q1 - Q2 + x(j);
    Q2 = Q1;
    Q1 = Q0;
end
A = 2*sqrt(Q1^2 + Q2^2 - Q1 * Q2 * coeff)/N;
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

              GOERTZELGH1TEST

goertzelgh1test uses goertzelgh1 to test for
the prescence of a single freqency in a set of
noisy signals.

% function [f1,A] = goertzelgh1test(N,Fs,SNR0,f0,ifsf,dthrsh)
%
% N Input signal lengths
% Fs Signal time sampling frequency
% SNR0 Signal-to-Noise ratio
%
% f0 Signal frequency to be detected
% ifsf Increased frequency sampling factor
% dthrsh Signal detection threshold (0 < dtrrsh <=1 )
%
% f1 Output frequencies within (approximately) the
% interval f0-Fs/N <= f1 <= f0+Fs/N
% A Amplitude of the frequency components in f

clear all, close all, clc

N = 44
Fs = 44.1e3
SNR0 = 32 % ~15dB
f0 = 4.5e3
ifsf = 10
dthrsh = 0.9

df = Fs/N % 1002.3 FFT resolution
f = df*(0:N-1)'; % FFT sampling
[fmin kmin] = max(f(f<f0)) % [4009.1 4]
fmax = f(kmin+1) % 5001.4

df0 = df/ifsf % 100.23 Goertzel sampling
M0 = 1+ceil((fmax-fmin)/(2*df0)) % 6
f1 = f0 + df0*[-M0:M0]'; % Goertzel sampling

M = 2*M0+1 % 13
minmaxf1 = minmax(f1') % [3898.6 5101.4]

dt = 1/Fs % 2.2676e-005
t = dt*(0:N-1)';

A = zeros(M,1);
for k = 1:M

    %Noisy test signal
    x0 = cos(2*pi*(f1(k)*t+rand));
    S0 = mean(x0.^2); % Average signal power
    x = x0 + sqrt(S0/SNR0)*randn(N,1);
    % Detection result
    A(k,1) = goertzelgh1(x,Fs,f0);

end
Athrsh = dthrsh*max(A)
summary = [ f1 A A>=Athrsh]

figure
plot(f1,A,'.-')
hold on
plot(f1,dthrsh*ones(size(f1)),'k--')
plot(f1(A>=Athrsh),A(A>=Athrsh),'r*')
legend('Output','Threshold','Detections')
ylim([ 0 1.5])
xlabel('Input frequency (Hz)')
ylabel('Detection Results')
title(['Frequency Detection Tests for f0 = ',num2str(f0),' Hz'])

On Jan 1, 10:35 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jan 1, 6:20 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> > On Dec 31 2010, 4:25 am, "Mario Rossi" <non...@hotmail.com> wrote:
> > -----SNIP
> > > Goertzel:
> > > for example, I'd like to use Goertzel, in a for cicle, to
> > > calculate all the single spectrum components in a range from
> > > 0,01 Hz and 5 Hz.
> > > But if I pass to the algorithm the non integer values (all
> > > values apart from 1, 2, 3, 4 ,5 Hz!!!) it returns an error!
> > > It says that it accepts only integer values..
> > > Do you know a way to avoid this problem?
>
> > I just modified Maurice Givens' codes in
>
> >http://groups.google.com/group/comp.dsp/msg/14c9c255d1005396
>
> > Will post them after dinner.
>
>                   GOERTZELGH1
>
> goertzelgh1 is used to detect a single frequency:
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> function A = goertzelgh1(x,Fs,f0)
>
> % Modified code of Maurice Givens
> %http://groups.google.com/group/comp.dsp/msg/...
> % 14c9c255d1005396

My modifications yield the code given in wikipedia.

> % x  = input vector
> % Fs = input sampling frequency
> % f0 = single test frequency
> % A  = test frequency amplplitude
> %
> % Greg Heath 1 Jan 2010
>
> if f0 >= Fs/2
>     error( 'f0 >= Nyquist frequency')
> end
>
> N     = length(x);
> coeff = 2 * cos(2*pi*f0/Fs);
> Q1 = 0; Q2 = 0;
> for j = 1:N
>     Q0 = coeff * Q1 - Q2 + x(j);
>     Q2 = Q1;
>     Q1 = Q0;
> end
> A = 2*sqrt(Q1^2 + Q2^2 - Q1 * Q2 * coeff)/N;
> return

% If used in a loop for detection, it may be
% more efficient to construct
% P = Q1^2 + Q2^2 - Q1 * Q2 * coeff;
% and use N to normalize the threshold.

Hope this helps.

Greg

Greg Heath <heath@alumni.brown.edu> wrote in message <0c1a17cb-ffaf-4468-8d65-4915864f8516@s9g2000vby.googlegroups.com>...
> On Jan 1, 10:35 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> > On Jan 1, 6:20 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> > > On Dec 31 2010, 4:25 am, "Mario Rossi" <non...@hotmail.com> wrote:
> > > -----SNIP
> > > > Goertzel:
> > > > for example, I'd like to use Goertzel, in a for cicle, to
> > > > calculate all the single spectrum components in a range from
> > > > 0,01 Hz and 5 Hz.
> > > > But if I pass to the algorithm the non integer values (all
> > > > values apart from 1, 2, 3, 4 ,5 Hz!!!) it returns an error!
> > > > It says that it accepts only integer values..
> > > > Do you know a way to avoid this problem?
> >
> > > I just modified Maurice Givens' codes in
> >
> > >http://groups.google.com/group/comp.dsp/msg/14c9c255d1005396
> >
> > > Will post them after dinner.
> >
> >                   GOERTZELGH1
> >
> > goertzelgh1 is used to detect a single frequency:
> >
> > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> >
> > function A = goertzelgh1(x,Fs,f0)
> >
> > % Modified code of Maurice Givens
> > %http://groups.google.com/group/comp.dsp/msg/...
> > % 14c9c255d1005396
>
> My modifications yield the code given in wikipedia.
>
> > % x  = input vector
> > % Fs = input sampling frequency
> > % f0 = single test frequency
> > % A  = test frequency amplplitude
> > %
> > % Greg Heath 1 Jan 2010
> >
> > if f0 >= Fs/2
> >     error( 'f0 >= Nyquist frequency')
> > end
> >
> > N     = length(x);
> > coeff = 2 * cos(2*pi*f0/Fs);
> > Q1 = 0; Q2 = 0;
> > for j = 1:N
> >     Q0 = coeff * Q1 - Q2 + x(j);
> >     Q2 = Q1;
> >     Q1 = Q0;
> > end
> > A = 2*sqrt(Q1^2 + Q2^2 - Q1 * Q2 * coeff)/N;
> > return
>
> % If used in a loop for detection, it may be
> % more efficient to construct
> % P = Q1^2 + Q2^2 - Q1 * Q2 * coeff;
> % and use N to normalize the threshold.
>
> Hope this helps.
>
> Greg

Thank you very much,
I'll try with you suggestions in these days..

Bye and happy new year!!

On Jan 1, 10:35 pm, Greg Heath <he...@alumni.brown.edu> wrote:

> % A = test frequency amplplitude

A = test frequency amplitude

> plot(f1,dthrsh*ones(size(f1)),'k--')

plot(f1,Athrsh*ones(size(f1)),'k--')

Hope this helps.

Greg

On Jan 2, 3:02 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jan 1, 10:35 pm, GregHeath<he...@alumni.brown.edu> wrote:
>
>
>
>
>
> > On Jan 1, 6:20 pm, GregHeath<he...@alumni.brown.edu> wrote:
> > > On Dec 31 2010, 4:25 am, "Mario Rossi" <non...@hotmail.com> wrote:
> > > -----SNIP
> > > > Goertzel:
> > > > for example, I'd like to use Goertzel, in a for cicle, to
> > > > calculate all the single spectrum components in a range from
> > > > 0,01 Hz and 5 Hz.
> > > > But if I pass to the algorithm the non integer values (all
> > > > values apart from 1, 2, 3, 4 ,5 Hz!!!) it returns an error!
> > > > It says that it accepts only integer values..
> > > > Do you know a way to avoid this problem?
>
> > > I just modified Maurice Givens' codes in
>
> > >http://groups.google.com/group/comp.dsp/msg/14c9c255d1005396
>
> > > Will post them after dinner.
>
> >                   GOERTZELGH1
>
> > goertzelgh1 is used to detect a single frequency:
>
> > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> > function A = goertzelgh1(x,Fs,f0)
>
> > % Modified code of Maurice Givens
> > %http://groups.google.com/group/comp.dsp/msg/...
> > % 14c9c255d1005396
>
> My modifications yield the code given in wikipedia.
>
>
>
>
>
> > % x  = input vector
> > % Fs = input sampling frequency
> > % f0 = single test frequency
> > % A  = test frequency amplplitude
> > %
> > % GregHeath1 Jan 2010
>
> > if f0 >= Fs/2
> >     error( 'f0 >= Nyquist frequency')
> > end
>
> > N     = length(x);
> > coeff = 2 * cos(2*pi*f0/Fs);
> > Q1 = 0; Q2 = 0;
> > for j = 1:N
> >     Q0 = coeff * Q1 - Q2 + x(j);
> >     Q2 = Q1;
> >     Q1 = Q0;
> > end
> > A = 2*sqrt(Q1^2 + Q2^2 - Q1 * Q2 * coeff)/N;
> > return
>
> % If used in a loop for detection, it may be
> % more efficient to construct
> % P = Q1^2 + Q2^2 - Q1 * Q2 * coeff;
> % and use N to normalize the threshold.

absX = sqrt(P) % Fourier amplitude
PSD = P/(N*Fs) * Power Spectral density

if 0 < abs(f0) < Fs/2

  % Should combine f0 and -f0 contributions

  absX = 2*sqrt(P)
  PSD = 2*P/(N*Fs)

end


Hope this helps.

Greg

On Dec 31 2010, 4:25 am, "Mario Rossi" <non...@hotmail.com> wrote:
> Greg Heath <he...@alumni.brown.edu> wrote in message <94453d96-2561-4fdf-aad7-9bb451170...@c2g2000yqc.googlegroups.com>...
> > On Dec 9, 6:55 am, "Mario Rossi" <non...@hotmail.com> wrote:
> > > Hi,
>
> > > I want to analyze the first 5 Hz of a signal spectrum..I've used the fft to calculate the spectrum, but it is 70 Hz wide, so I have a lot of "useless" frequency..
>
> > > It's possible to "concentrate" the 1024 FFT points (number I'm actually using) in the first (about) 5 Hz?
> > > In this way I could get an higher frequency resolution, higher than the actual 0.1 Hz.
>
> > I'm confused.
>
> > To get a resolution of df = 0.1 Hz with N = 1024
> > points, you need a sampling frequency of
>
> > Fs = N*df = 102.4
> > dt = 1/Fs = 9.766 msec
> > tmax = (N-1)*dt = 9.99 sec
> > fmax = Fs/2 = 51.2 Hz
>
> > Where does "70 Hz" come from?
>
> > > I've tried theGoertzelalgorithm, but it works
> > > only with integer frequency indexes, so it doesn't
> > > work fine for me..

Use the algorithm given in Wikipedia which I
posted previously.


> About your second post..
> 70 Hz come from the length of the signal = 14 s, and
the number of samples (1000) in this signal.
> fs = 1/(14/1000) = 71 Hz (more or less)
> I was wrong when I said "actual resolution of 0.1 Hz",
> it is about 0.07 Hz. However the problem is still
> the same, I have to improve the resolution.

Do you mean resolution or frequency sampling interval?

> Goertzel:
> for example, I'd like to use Goertzel, in a for cicle,
to calculate all the single spectrum components in a
range from 0,01 Hz and 5 Hz.
> But if I pass to the algorithm the non integer values
(all values apart from 1, 2, 3, 4 ,5 Hz!!!) it returns
an error! It says that it accepts only integer values..
> Do you know a way to avoid this problem?

See above.

Are the following problem specifications correct?

% GIVEN:

1. x(t), 0 <= t <= tmax
1. dt = 1/Fs fixed by hardware
2. tmax = 14
3. N = 1000 % Not 1024 ?

% resulting in

dt = tmax/(N-1) % 0.014014
Fs = 1/dt % 71.357
df0 = Fs/N % 0.071357
T = N*dt % 14.014

% df0 = 1/T = 1/((tmax+dt)

% FIND

1. Frequency Spectrum
2. fmin = 0, fmax = 5
3. df = 0.01

% Questions

1. Is x real valued?
2. Is x a single frequency sinusoid with
   additive Gaussian noise?
3. Do you want PSD or Ampl & Phase?

% SOLUTIONS

1. To get true 0.01 resolution, either
   a. increase tmax and N @ Fs=constant
   b. Find MATLAB superresolution functions:
      lookfor music
2. To get 0.01 interpolation consider
   a. Zeropad and fft
   b. dft
   c. Goertzel
   d. Chirp z-transform
3. Corresponding calculation sizes (real x):
   a. zeropad and fft (df = Fs/Mz <= 0.01)
      Mz = 2^ceil(log2(Fs/df)) % 8192
   b. dft & goertzel
      Md = 1 + ceil((fmax-fmin)/df) % 501
   c. Chirp-z
      Same as goertzel ???

The important point here is that

Md = 501 >> ceil(log2(Mz)) = 13

Therefore, it should be faster to use
the zero-padded fft.

Hope this helps.

Greg

> See above.
>
> Are the following problem specifications correct?
>
> % GIVEN:
>
> 1. x(t), 0 <= t <= tmax
> 1. dt = 1/Fs fixed by hardware
> 2. tmax = 14
> 3. N = 1000 % Not 1024 ?
>
> % resulting in
>
> dt = tmax/(N-1) % 0.014014
> Fs = 1/dt % 71.357
> df0 = Fs/N % 0.071357
> T = N*dt % 14.014
>
> % df0 = 1/T = 1/((tmax+dt)
>
> % FIND
>
> 1. Frequency Spectrum
> 2. fmin = 0, fmax = 5
> 3. df = 0.01
>
> % Questions
>
> 1. Is x real valued?
> 2. Is x a single frequency sinusoid with
> additive Gaussian noise?
> 3. Do you want PSD or Ampl & Phase?
>
> % SOLUTIONS
>
> 1. To get true 0.01 resolution, either
> a. increase tmax and N @ Fs=constant
> b. Find MATLAB superresolution functions:
> lookfor music
> 2. To get 0.01 interpolation consider
> a. Zeropad and fft
> b. dft
> c. Goertzel
> d. Chirp z-transform
> 3. Corresponding calculation sizes (real x):
> a. zeropad and fft (df = Fs/Mz <= 0.01)
> Mz = 2^ceil(log2(Fs/df)) % 8192
> b. dft & goertzel
> Md = 1 + ceil((fmax-fmin)/df) % 501
> c. Chirp-z
> Same as goertzel ???
>
> The important point here is that
>
> Md = 501 >> ceil(log2(Mz)) = 13
>
> Therefore, it should be faster to use
> the zero-padded fft.
>
> Hope this helps.
>
> Greg

Thank you Greg, but why do you reply to my posts "randomly"? :D

However your previous version of Goertzel doesn't work fine for me, maybe I make something wrong, I don't know..

About the zero-padding... zero pad = interpolation... I'd like to use a method that is more accurate... after all, interpolation is an approximation...

However... I'm still working on Goertzel, then I'll post the results..
(I've already looked for music algorithm, but It's to advanced for me, I'm not able to use it :( )

Bye

On Jan 9, 6:49 am, "Mario Rossi" <non...@hotmail.com> wrote:
> > See above.
>
> > Are the following problem specifications correct?

PLEASE ANSWER

> > % GIVEN:
>
> > 1. x(t), 0 <= t <= tmax
> > 1. dt = 1/Fs fixed by hardware
> > 2. tmax = 14
> > 3. N    = 1000       % Not 1024 ?
>
> > % resulting in
>
> > dt  = tmax/(N-1)   % 0.014014
> > Fs  = 1/dt         % 71.357
> > df0 = Fs/N         % 0.071357
> > T   = N*dt         % 14.014
>
> > % df0 = 1/T = 1/((tmax+dt)
>
> > % FIND
>
> > 1. Frequency Spectrum
> > 2. fmin = 0, fmax = 5
> > 3. df = 0.01
>
> > % Questions
>
> > 1. Is x real valued?
> > 2. Is x a single frequency sinusoid with
> >    additive Gaussian noise?
> > 3. Do you want PSD or Ampl & Phase?

PLEASE ANSWER

> > % SOLUTIONS
>
> > 1. To get true 0.01 resolution, either
> >    a. increase tmax and N @ Fs=constant
> >    b. Find MATLAB superresolution functions:
> >       lookfor music
> > 2. To get 0.01 interpolation consider
> >    a. Zeropad and fft
> >    b. dft
> >    c. Goertzel
> >    d. Chirp z-transform
> > 3. Corresponding calculation sizes (real x):
> >    a. zeropad and fft (df = Fs/Mz <= 0.01)
> >       Mz = 2^ceil(log2(Fs/df))  % 8192
> >    b. dft & goertzel
> >       Md = 1 + ceil((fmax-fmin)/df) % 501
> >    c. Chirp-z
> >       Same as goertzel???
>
> > The important point here is that
>
> > Md = 501 >> ceil(log2(Mz)) = 13
>
> > Therefore, it should be faster to use
> > the zero-padded fft.
>
> > Hope this helps.
>
> > Greg
>
> Thank you Greg, but why do you reply to my posts "randomly"? :D

I use the Google Groups reader which has
the option of listing the posts either chronologically or in a tree
where each
post is followed by it's direct replies.

Unfortunately my short-term memory is decaying
with age. Consequently, I frequently review long threads in tree
format to keep from losing
track. As I review a post and think of something
new, I reply to that post even though it is not
the most recent.

> However your previous version of Goertzeldoesn't work fine for me, maybe I make something wrong, I >don't know..

Why do you think it is not working as it should?
What data and parameters are you using to test it?
Did you try sinusoid + noise?

> About the zero-padding... zero pad = interpolation... I'd like to use a method that is more accurate... after all, interpolation is an >approximation...

Since zero-padding is the most common use of
the fft, it should be instructive to compare
zero padding with dft and goertzel. By the way,
if the spectra are evenly spaced, there is
a faster way to calculate the dft (Horner's method)
I will dig up some timing tests for fft vs Horner and matrix dft.

Again, is my above specification of the problem
correct?

> However... I'm still working on Goertzel, then >I'll post the results..

Again, using smaller frequency spacing in dft and Goertzel and does
not improve resolution.

> (I've already looked for music algorithm, but It's >to advanced for me, I'm not able to use it :(  )

I've never used the superresolution algorithms.
Right now I want to look at zoom-fft, pruned-fft
and chirp-z.

Are you interested in amplitude and phase
(GoertzelAgh) or is PSD (GoertzelPgh) sufficient?

Hope this helps.

Greg

On Jan 9, 12:22 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> On Jan 9, 6:49 am, "Mario Rossi" <non...@hotmail.com> wrote:
>
> > > See above.
>
> > > Are the following problem specifications correct?
>
> PLEASE ANSWER
>
> > > % GIVEN:
>
> > > 1. x(t), 0 <= t <= tmax
> > > 1. dt = 1/Fs fixed by hardware
> > > 2. tmax = 14
> > > 3. N    = 1000       % Not 1024 ?
>
> > > % resulting in
>
> > > dt  = tmax/(N-1)   % 0.014014
> > > Fs  = 1/dt         % 71.357
> > > df0 = Fs/N         % 0.071357
> > > T   = N*dt         % 14.014
>
> > > % df0 = 1/T = 1/((tmax+dt)
>
> > > % FIND
>
> > > 1. Frequency Spectrum
> > > 2. fmin = 0, fmax = 5
> > > 3. df = 0.01
>
> > > % Questions
>
> > > 1. Is x real valued?
> > > 2. Is x a single frequency sinusoid with
> > >    additive Gaussian noise?
> > > 3. Do you want PSD or Ampl & Phase?
>
> PLEASE ANSWER
>
> > > % SOLUTIONS
>
> > > 1. To get true 0.01 resolution, either
> > >    a. increase tmax and N @ Fs=constant
> > >    b. Find MATLAB superresolution functions:
> > >       lookfor music
> > > 2. To get 0.01 interpolation consider
> > >    a. Zeropad and fft
> > >    b. dft
> > >    c. Goertzel
> > >    d. Chirp z-transform
> > > 3. Corresponding calculation sizes (real x):
> > >    a. zeropad and fft (df = Fs/Mz <= 0.01)
> > >       Mz = 2^ceil(log2(Fs/df))  % 8192
> > >    b. dft & goertzel
> > >       Md = 1 + ceil((fmax-fmin)/df) % 501
> > >    c. Chirp-z
> > >       Same as goertzel???
>
> > > The important point here is that
>
> > > Md = 501 >> ceil(log2(Mz)) = 13
>
> > > Therefore, it should be faster to use
> > > the zero-padded fft.

I take that back if you are willing to use
C-Mex.

> > Thank you Greg, but why do you reply to my posts
>>"randomly"? :D
>
> I use the Google Groups reader which has
> the option of listing the posts either
> chronologically or in a tree where each
> post is followed by it's direct replies.
>
> Unfortunately my short-term memory is decaying
> with age. Consequently, I frequently review
> long threads in tree format to keep from losing
> track. As I review a post and think of something
> new, I reply to that post even though it is not
> the most recent.
>
> > However your previous version of Goertzel doesn't
> >work fine for me, maybe I make something wrong, I
>>don't know..
>
> Why do you think it is not working as it should?
> What data and parameters are you using to test it?
> Did you try sinusoid + noise?
>
> > About the zero-padding... zero pad = interpolation
... I'd like to use a method that is more accurate...
>>after all, interpolation is an >approximation...
>
> Since zero-padding is the most common use of
> the fft, it should be instructive to compare
> zero padding with dft and goertzel. By the way,
> if the spectra are evenly spaced, there is
> a faster way to calculate the dft (Horner's method)
> I will dig up some timing tests for fft vs Horner
> and matrix dft.
>
> Again, is my above specification of the problem
> correct?
>
> > However... I'm still working on Goertzel, then
>I'll post the results..
>
> Again, using smaller frequency spacing in dft and
> Goertzel and does not improve resolution.
>
> > (I've already looked for music algorithm, but
> > It's >to advanced for me, I'm not able to use it :(  )
>
> I've never used the superresolution algorithms.
> Right now I want to look at zoom-fft, pruned-fft
> and chirp-z.

So far:

Zoom FFT:
1. Multiply the data by a sinusoid with the midband
   frequency f0. This yields a difference frequency
   spectrum component centered at DC and a sum frequency
  spectrum component centered at 2*f0
2. LPF to get rid of the sum frequencies and higher
   difference frequencies
3. Sample the result and fft
4. Ordinarily the resampling is done at a lower
   Fs appropriate to the new maximum frequency.
   However, you said your available Fs was not variable.
5. Since you are interested in baseband, you can
   skip step1.

This corresponds to the replies of Nic and Matt.

Pruned FFT
1. Take fft
2. Remove desired frequencies
3. Reconstruct with idft (not ifft) or inverse
   Goetzel.
4. Use Mexed C functions to speed up the
   inversions in step 3.

So far it looks like the zoom-fft makes the most
sense.

Again:

Are you interested in amplitude and phase
(e.g., GoertzelAgh) or is PSD (e.g., GoertzelPgh)
sufficient?

Hope this helps.

Greg

Greg Heath <heath@alumni.brown.edu> wrote in message <ba37cb61-22aa-4c53-a593-e0036d75dcfe@i41g2000vbn.googlegroups.com>...
> On Jan 9, 12:22 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> > On Jan 9, 6:49 am, "Mario Rossi" <non...@hotmail.com> wrote:
> >
> > > > See above.
> >
> > > > Are the following problem specifications correct?
> >
> > PLEASE ANSWER
> >
> > > > % GIVEN:
> >
> > > > 1. x(t), 0 <= t <= tmax
> > > > 1. dt = 1/Fs fixed by hardware
> > > > 2. tmax = 14
> > > > 3. N    = 1000       % Not 1024 ?
> >
> > > > % resulting in
> >
> > > > dt  = tmax/(N-1)   % 0.014014
> > > > Fs  = 1/dt         % 71.357
> > > > df0 = Fs/N         % 0.071357
> > > > T   = N*dt         % 14.014
> >
> > > > % df0 = 1/T = 1/((tmax+dt)
> >
> > > > % FIND
> >
> > > > 1. Frequency Spectrum
> > > > 2. fmin = 0, fmax = 5
> > > > 3. df = 0.01
> >
> > > > % Questions
> >
> > > > 1. Is x real valued?
> > > > 2. Is x a single frequency sinusoid with
> > > >    additive Gaussian noise?
> > > > 3. Do you want PSD or Ampl & Phase?
> >
> > PLEASE ANSWER
> >
> > > > % SOLUTIONS
> >
> > > > 1. To get true 0.01 resolution, either
> > > >    a. increase tmax and N @ Fs=constant
> > > >    b. Find MATLAB superresolution functions:
> > > >       lookfor music
> > > > 2. To get 0.01 interpolation consider
> > > >    a. Zeropad and fft
> > > >    b. dft
> > > >    c. Goertzel
> > > >    d. Chirp z-transform
> > > > 3. Corresponding calculation sizes (real x):
> > > >    a. zeropad and fft (df = Fs/Mz <= 0.01)
> > > >       Mz = 2^ceil(log2(Fs/df))  % 8192
> > > >    b. dft & goertzel
> > > >       Md = 1 + ceil((fmax-fmin)/df) % 501
> > > >    c. Chirp-z
> > > >       Same as goertzel???
> >
> > > > The important point here is that
> >
> > > > Md = 501 >> ceil(log2(Mz)) = 13
> >
> > > > Therefore, it should be faster to use
> > > > the zero-padded fft.
>
> I take that back if you are willing to use
> C-Mex.
>
> > > Thank you Greg, but why do you reply to my posts
> >>"randomly"? :D
> >
> > I use the Google Groups reader which has
> > the option of listing the posts either
> > chronologically or in a tree where each
> > post is followed by it's direct replies.
> >
> > Unfortunately my short-term memory is decaying
> > with age. Consequently, I frequently review
> > long threads in tree format to keep from losing
> > track. As I review a post and think of something
> > new, I reply to that post even though it is not
> > the most recent.
> >
> > > However your previous version of Goertzel doesn't
> > >work fine for me, maybe I make something wrong, I
> >>don't know..
> >
> > Why do you think it is not working as it should?
> > What data and parameters are you using to test it?
> > Did you try sinusoid + noise?
> >
> > > About the zero-padding... zero pad = interpolation
> ... I'd like to use a method that is more accurate...
> >>after all, interpolation is an >approximation...
> >
> > Since zero-padding is the most common use of
> > the fft, it should be instructive to compare
> > zero padding with dft and goertzel. By the way,
> > if the spectra are evenly spaced, there is
> > a faster way to calculate the dft (Horner's method)
> > I will dig up some timing tests for fft vs Horner
> > and matrix dft.
> >
> > Again, is my above specification of the problem
> > correct?
> >
> > > However... I'm still working on Goertzel, then
> >I'll post the results..
> >
> > Again, using smaller frequency spacing in dft and
> > Goertzel and does not improve resolution.
> >
> > > (I've already looked for music algorithm, but
> > > It's >to advanced for me, I'm not able to use it :(  )
> >
> > I've never used the superresolution algorithms.
> > Right now I want to look at zoom-fft, pruned-fft
> > and chirp-z.
>
> So far:
>
> Zoom FFT:
> 1. Multiply the data by a sinusoid with the midband
> frequency f0. This yields a difference frequency
> spectrum component centered at DC and a sum frequency
> spectrum component centered at 2*f0
> 2. LPF to get rid of the sum frequencies and higher
> difference frequencies
> 3. Sample the result and fft
> 4. Ordinarily the resampling is done at a lower
> Fs appropriate to the new maximum frequency.
> However, you said your available Fs was not variable.
> 5. Since you are interested in baseband, you can
> skip step1.
>
> This corresponds to the replies of Nic and Matt.
>
> Pruned FFT
> 1. Take fft
> 2. Remove desired frequencies
> 3. Reconstruct with idft (not ifft) or inverse
> Goetzel.
> 4. Use Mexed C functions to speed up the
> inversions in step 3.
>
> So far it looks like the zoom-fft makes the most
> sense.
>
> Again:
>
> Are you interested in amplitude and phase
> (e.g., GoertzelAgh) or is PSD (e.g., GoertzelPgh)
> sufficient?
>
> Hope this helps.
>
> Greg

Hi Greg,

I'm in a hurry so I'll reply to you more accurately in these days..
However:
- x is real.
- it hasn't got only one sinusoid and gaussian noise, but several frequencies+noise (it's a signal received from an antenna), my goal is to find if it's present or not a frequency component at a certain frequency.
- I don't need phase information, I need only amplitude.

( sorry for my bad english :) )

Thanks,

Mario

On Jan 10, 4:05 am, "Mario Rossi" <non...@hotmail.com> wrote:
> Greg Heath <he...@alumni.brown.edu> wrote in message <ba37cb61-22aa-4c53-a593-e0036d75d...@i41g2000vbn.googlegroups.com>...
> > On Jan 9, 12:22 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> > > On Jan 9, 6:49 am, "Mario Rossi" <non...@hotmail.com> wrote:
>
> > > > > See above.
>
> > > > > Are the following problem specifications correct?
>
> > > PLEASE ANSWER
>
> > > > > % GIVEN:
>
> > > > > 1. x(t), 0 <= t <= tmax
> > > > > 1. dt = 1/Fs fixed by hardware
> > > > > 2. tmax = 14
> > > > > 3. N = 1000 % Not 1024 ?
>
> > > > > % resulting in
>
> > > > > dt = tmax/(N-1) % 0.014014
> > > > > Fs = 1/dt % 71.357
> > > > > df0 = Fs/N % 0.071357
> > > > > T = N*dt % 14.014
>
> > > > > % df0 = 1/T = 1/((tmax+dt)
>
> > > > > % FIND
>
> > > > > 1. Frequency Spectrum
> > > > > 2. fmin = 0, fmax = 5
> > > > > 3. df = 0.01
>
> > > > > % Questions
>
> > > > > 1. Is x real valued?
> > > > > 2. Is x a single frequency sinusoid with
> > > > > additive Gaussian noise?
> > > > > 3. Do you want PSD or Ampl & Phase?
>
> > > PLEASE ANSWER
>
> > > > > % SOLUTIONS
>
> > > > > 1. To get true 0.01 resolution, either
> > > > > a. increase tmax and N @ Fs=constant
> > > > > b. Find MATLAB superresolution functions:
> > > > > lookfor music
> > > > > 2. To get 0.01 interpolation consider
> > > > > a. Zeropad and fft
> > > > > b. dft
> > > > > c. Goertzel
> > > > > d. Chirp z-transform
> > > > > 3. Corresponding calculation sizes (real x):
> > > > > a. zeropad and fft (df = Fs/Mz <= 0.01)
> > > > > Mz = 2^ceil(log2(Fs/df)) % 8192
> > > > > b. dft & goertzel
> > > > > Md = 1 + ceil((fmax-fmin)/df) % 501
> > > > > c. Chirp-z
> > > > > Same as goertzel???
>
> > > > > The important point here is that
>
> > > > > Md = 501 >> ceil(log2(Mz)) = 13
>
> > > > > Therefore, it should be faster to use
> > > > > the zero-padded fft.
>
> > I take that back if you are willing to use
> > C-Mex.
>
> > > > Thank you Greg, but why do you reply to my posts
> > >>"randomly"? :D
>
> > > I use the Google Groups reader which has
> > > the option of listing the posts either
> > > chronologically or in a tree where each
> > > post is followed by it's direct replies.
>
> > > Unfortunately my short-term memory is decaying
> > > with age. Consequently, I frequently review
> > > long threads in tree format to keep from losing
> > > track. As I review a post and think of something
> > > new, I reply to that post even though it is not
> > > the most recent.
>
> > > > However your previous version of Goertzel doesn't
> > > >work fine for me, maybe I make something wrong, I
> > >>don't know..
>
> > > Why do you think it is not working as it should?
> > > What data and parameters are you using to test it?
> > > Did you try sinusoid + noise?
>
> > > > About the zero-padding... zero pad = interpolation
> > ... I'd like to use a method that is more accurate...
> > >>after all, interpolation is an >approximation...
>
> > > Since zero-padding is the most common use of
> > > the fft, it should be instructive to compare
> > > zero padding with dft and goertzel. By the way,
> > > if the spectra are evenly spaced, there is
> > > a faster way to calculate the dft (Horner's method)
> > > I will dig up some timing tests for fft vs Horner
> > > and matrix dft.
>
> > > Again, is my above specification of the problem
> > > correct?
>
> > > > However... I'm still working on Goertzel, then
> > >I'll post the results..
>
> > > Again, using smaller frequency spacing in dft and
> > > Goertzel and does not improve resolution.
>
> > > > (I've already looked for music algorithm, but
> > > > It's >to advanced for me, I'm not able to use it :( )
>
> > > I've never used the superresolution algorithms.
> > > Right now I want to look at zoom-fft, pruned-fft
> > > and chirp-z.
>
> > So far:
>
> > Zoom FFT:
> > 1. Multiply the data by a sinusoid with the midband
> > frequency f0. This yields a difference frequency
> > spectrum component centered at DC and a sum frequency
> > spectrum component centered at 2*f0
> > 2. LPF to get rid of the sum frequencies and higher
> > difference frequencies
> > 3. Sample the result and fft
> > 4. Ordinarily the resampling is done at a lower
> > Fs appropriate to the new maximum frequency.
> > However, you said your available Fs was not variable.
> > 5. Since you are interested in baseband, you can
> > skip step1.
>
> > This corresponds to the replies of Nic and Matt.
>
> > Pruned FFT
> > 1. Take fft
> > 2. Remove desired frequencies
> > 3. Reconstruct with idft (not ifft) or inverse
> > Goetzel.
> > 4. Use Mexed C functions to speed up the
> > inversions in step 3.
>
> > So far it looks like the zoom-fft makes the most
> > sense.
>
> > Again:
>
> > Are you interested in amplitude and phase
> > (e.g., GoertzelAgh) or is PSD (e.g., GoertzelPgh)
> > sufficient?
>
> > Hope this helps.
>
> > Greg
>
> Hi Greg,
>
> I'm in a hurry so I'll reply to you more accurately in these days..
> However:
> - x is real.
> - it hasn't got only one sinusoid and gaussian noise, but several frequencies+noise (it's a signal received from an antenna), my goal is to find if it's present or not a frequency component at a certain frequency.
> - I don't need phase information, I need only amplitude.
>
> ( sorry for my bad english :) )
>
> Thanks,
>
> Mario- Hide quoted text -
>
> - Show quoted text -

How many sinusoids in the 0-5 Hz range?
How much can their frequencies differ?
How much can their magnitudes differ?
What is the range of average noise power?

Hope this helps.

Greg

> How many sinusoids in the 0-5 Hz range?
> How much can their frequencies differ?
> How much can their magnitudes differ?
> What is the range of average noise power?
>
> Hope this helps.
>
> Greg

Hi greg,

I've spoken with my collegues and we have decided to quit this approach to the problem... I've seen that the higher resolution won't give an improvement to this project..

However, just to answer to your questions:
- I have a received signal that is the result of the interation of multiple sent pulses with the environment.
- the noise is (approximately) white.

Thank you for your useful advices.

Bye

On Jan 11, 3:26 am, "Mario Rossi" <non...@hotmail.com> wrote:
> > How many sinusoids in the 0-5 Hz range?
> > How much can their frequencies differ?
> > How much can their magnitudes differ?
> > What is the range of average noise power?
>
> > Hope this helps.
>
> > Greg
>
> Hi greg,
>
> I've spoken with my collegues and we have decided to quit this approach to the problem... I've seen that the higher resolution won't give an improvement to this project..
>
> However, just to answer to your questions:
> - I have a received signal that is the result of the interation of multiple sent pulses with the environment.
> - the noise is (approximately) white.
>
> Thank you for your useful advices.

You are welcome.

What can uou say about the following questions:

 How many sinusoids in the 0-5 Hz range?
 How much can their frequencies differ?
 How much can their magnitudes differ?
 What is the range of average noise power?

TIA

Greg

Mario Rossi wrote:
>> How many sinusoids in the 0-5 Hz range?
>> How much can their frequencies differ?
>> How much can their magnitudes differ?
>> What is the range of average noise power?
>>
>> Hope this helps.
>>
>> Greg
>
> Hi greg,
> I've spoken with my collegues and we have decided to quit this approach
> to the problem... I've seen that the higher resolution won't give an
> improvement to this project..
> However, just to answer to your questions:
> - I have a received signal that is the result of the interation of
> multiple sent pulses with the environment.
> - the noise is (approximately) white.
...

Besides the queries from Greg, white noise will be amenable to averaging...

--

> How many sinusoids in the 0-5 Hz range?
I don't know, the signal I'm analysing is the result of a lot of components.

> How much can their frequencies differ?
> How much can their magnitudes differ?
> What is the range of average noise power?
the average noise is present at all the frequencies...


> TIA
>
> Greg

Don't waste further time for me Greg, I'm not going to use this approach to the problem anymore :)

Thanks,

Mario

On Jan 13, 7:32 am, "Mario Rossi" <non...@hotmail.com> wrote:
> >  How many sinusoids in the 0-5 Hz range?
>
> I don't know, the signal I'm analysing is the result of a lot of components.
>
> >  How much can their frequencies differ?
> >  How much can their magnitudes differ?
> >  What is the range of average noise power?
>
> the average noise is present at all the frequencies...
>
> > TIA
>
> > Greg
>
> Don't waste further time for me Greg, I'm not going to >use this approach to the problem anymore :)

I don't consider the time wasted.

I need interesting problems to keep my
brain from rotting (age and health have
put a huge damper on sports and carousing).

Using RW data and/or parameters feels
better than arbitraily choosing some
numbers out of mid-air.

Greg

> I don't consider the time wasted.
>
> I need interesting problems to keep my
> brain from rotting (age and health have
> put a huge damper on sports and carousing).
>
> Using RW data and/or parameters feels
> better than arbitraily choosing some
> numbers out of mid-air.
>
> Greg

Ok.
However I don't have specific data.
I receive data to process from a radar antenna, and so se signal is a mess..
So I haven't got any information about frequency...

Mario

On Jan 14, 3:22 am, "Mario Rossi" <non...@hotmail.com> wrote:
> > I don't consider the time wasted.
>
> > I need interesting problems to keep my
> > brain from rotting (age and health have
> > put a huge damper on sports and carousing).
>
> > Using RW data and/or parameters feels
> > better than arbitraily choosing some
> > numbers out of mid-air.
>
> > Greg
>
> Ok.
> However I don't have specific data.
> I receive data to process from a radar antenna, and so se signal is a mess..
> So I haven't got any information about frequency...
>
> Mario

OK. Currently using Fs = 70, N = 1000, df = Fs/N = 0.07,
f0 = 2.5 in three different scenarios with three frequencies each:

1. f0-df/2, f0, f0+df/2
2. f0-df, f0, f0+df
3. f0-3*df/2, f0, f0+3*df/2

Greg