Re: Real-data DFT in Octave

[Sergei Steshenko]

----- Original Message -----
> From: Torbjörn Rathsman <address@hidden>
> To: Sergei Steshenko <address@hidden>; "address@hidden" <address@hidden>
> Cc: 
> Sent: Saturday, October 6, 2012 9:05 PM
> Subject: Re: Real-data DFT in Octave
> 
> Sergei Steshenko skrev 2012-10-06 20:47:
>> 
>> 
>> 
>>  ----- Original Message -----
>>>  From: Torbjörn Rathsman <address@hidden>
>>>  To: Sergei Steshenko <address@hidden>; 
> "address@hidden" <address@hidden>
>>>  Cc:
>>>  Sent: Saturday, October 6, 2012 8:25 PM
>>>  Subject: Re: Real-data DFT in Octave
>>> 
>>>  Sergei Steshenko skrev 2012-10-06 20:06:
>>>> 
>>>> 
>>>>    ----- Original Message -----
>>>>>    From: Torbjörn Rathsman <address@hidden>
>>>>>    To: Sergei Steshenko <address@hidden>
>>>>>    Cc:
>>>>>    Sent: Saturday, October 6, 2012 7:36 PM
>>>>>    Subject: Re: Real-data DFT in Octave
>>>>> 
>>>>    [snip]
>>>>>    I get imaginary part that is in the same order of magnitude 
> as the real
>>>>>    part. Note: my N is odd, but that should only result in a 
> discarded
>>>>>    sample or?
>>>>    For odd N you do the following:
>>>> 
>>>>    1) leave (N+1)/2 frequency bin intact;
>>>>    2) take the first (N-1)/2 frequency bins and mirror them around 
> (N+1)/2
>>>  frequency bin.
>>>>    In such a manner nothing is lost/discarded, everything is 
> mathematically
>>>  strictly invertible.
>>>>    E.g. for N=5 you have in Octave terms:
>>>> 
>>>>    1 2 3 2' 1'
>>>> 
>>>>    where "'" means complex conjugate; since 
> '1' is DC
>>>  1' == 1, but who cares, i.e. it's easier to write mirroring not 
> taking
>>>  this into account - complex conjugate of a real component is still the 
> same
>>>  real.
>>>>    Regards,
>>>>       Sergei.
>>>> 
>>>> 
>>>  Even N also gives me wrong result.
>>> 
>>>       Y=fft(y);
>>>       length(y) %31047
>>>       Y=Y(1:floor(length(Y)/2));
>>>       length(Y) %15523
>>>       for k=1:length(Y)
>>>           if(k<=1024)
>>>               Y(k)=(k-1)*(k-1)*Y(k)/N; %diffrentiate twice
>>>           else
>>>               Y(k)=0;
>>>           end
>>>       end
>>>       Y=[Y; conj(flipud(Y))];
>>>       y=ifft(y);
>>>       plot(real(y),imag(y));
>> 
>>  What are you trying to achieve ? If you have spectrum, why won't use it 
> to calculate derivative ?
>> 
>> 
>>  Regards,
>>     Sergei.
>> 
>> 
> I am trying to apply a filter that is k*k for k<1024 which is somewhat 
> equivalent to diffrentiate the input twice for low frequencies. Also I 
> want to get rid of high frequency noise. However, the filter does not work
> 
> Now y is of even length:
> 
>      Y=fft(y);
>      length(y)
>      Y=Y(1:floor(length(Y)/2));
>      length(Y)
> 
>      Y=[Y; conj(flipud(Y))];
>      y=ifft(Y); %I missed capital Y here but
>      plot(real(y),imag(y)); %still huge imaginary part.

Sorry, I screwed up with off by one error.

Here is my screen session:

"
octave:1> signal = [1 2 3 4 5 6 7 8]
signal =

   1   2   3   4   5   6   7   8

octave:2> spectrum = fft(signal)
spectrum =

   36.0000 +  0.0000i   -4.0000 +  9.6569i   -4.0000 +  4.0000i   -4.0000 +  
1.6569i   -4.0000 +  0.0000i   -4.0000 -  1.6569i   -4.0000 -  4.0000i   
-4.0000 -  9.6569i

octave:3> spectrum(end/2+2:end) = conj(fliplr(spectrum(2:end/2)))
spectrum =

   36.0000 +  0.0000i   -4.0000 +  9.6569i   -4.0000 +  4.0000i   -4.0000 +  
1.6569i   -4.0000 +  0.0000i   -4.0000 -  1.6569i   -4.0000 -  4.0000i   
-4.0000 -  9.6569i

octave:4> ifft(spectrum)
ans =

   1   2   3   4   5   6   7   8

octave:5>
"

- please pay attention to mirroring:

spectrum(end/2+2:end) = conj(fliplr(spectrum(2:end/2)))
.

I.e. the DC component should not take part in the mirroring. And in this case 
of even N mirroring takes part around N/2+1 bin which is left intact.

Regards,
  Sergei.






>