V_IRFFT Inverse fft of a conjugate symmetric spectrum X=(Y,N,D) Inputs: Y(M) The first half of a complex spectrum N The number of output points to generate (default: 2M-2) D The dimension along which to perorm the transform (default: first non-singleton dimension of Y) Outputs: X(N) Real inverse dft of Y This routine calculates the inverse DFT of a conjugate-symmetric to give a real-valued output of dimension N. Only the first half of the spectrum need be supplied: if N is even, this includes the Nyquist term and is of dimension M=N/2 + 1 whereas if N is odd then there is no Nyquist term and the input is of dimension M=(N+1)/2. The DC and (if present) Nyquist values are forced to be real before performing the transform. Note that the default value of N is always even so that N must be given explicitly if it is odd. See also the forward transform: RFFT
0001 function x=v_irfft(y,n,d) 0002 %V_IRFFT Inverse fft of a conjugate symmetric spectrum X=(Y,N,D) 0003 % 0004 % Inputs: Y(M) The first half of a complex spectrum 0005 % N The number of output points to generate (default: 2M-2) 0006 % D The dimension along which to perorm the transform 0007 % (default: first non-singleton dimension of Y) 0008 % 0009 % Outputs: X(N) Real inverse dft of Y 0010 % 0011 % This routine calculates the inverse DFT of a conjugate-symmetric to give a real-valued 0012 % output of dimension N. Only the first half of the spectrum need be supplied: if N is even, 0013 % this includes the Nyquist term and is of dimension M=N/2 + 1 whereas if N is odd then there is 0014 % no Nyquist term and the input is of dimension M=(N+1)/2. 0015 % The DC and (if present) Nyquist values are forced to be real before performing the transform. 0016 % Note that the default value of N is always even so that N must be given explicitly 0017 % if it is odd. 0018 % 0019 % See also the forward transform: RFFT 0020 0021 % Copyright (C) Mike Brookes 2009 0022 % Version: $Id: v_irfft.m 10865 2018-09-21 17:22:45Z dmb $ 0023 % 0024 % VOICEBOX is a MATLAB toolbox for speech processing. 0025 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0026 % 0027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0028 % This program is free software; you can redistribute it and/or modify 0029 % it under the terms of the GNU General Public License as published by 0030 % the Free Software Foundation; either version 2 of the License, or 0031 % (at your option) any later version. 0032 % 0033 % This program is distributed in the hope that it will be useful, 0034 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0035 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0036 % GNU General Public License for more details. 0037 % 0038 % You can obtain a copy of the GNU General Public License from 0039 % http://www.gnu.org/copyleft/gpl.html or by writing to 0040 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0041 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0042 0043 s=size(y); 0044 ps=prod(s); 0045 ns=length(s); 0046 if ps==1 0047 x=y; 0048 else 0049 if nargin <3 || isempty(d) 0050 d=find(s>1,1); 0051 end 0052 m=s(d); 0053 k=ps/m; % number of fft's to do 0054 if d==1 0055 v=reshape(y,m,k); 0056 else 0057 v=reshape(permute(y,[d:ns 1:d-1]),m,k); 0058 end 0059 if nargin<2 || isempty(n) 0060 n=2*m-2; % default output length 0061 else 0062 mm=1+fix(n/2); % expected input length 0063 if mm>m v=[v; zeros(mm-m,k)]; % zero pad 0064 elseif mm<m v(mm+1:m,:)=[]; % or truncate 0065 end 0066 m=mm; 0067 end 0068 v(1,:)=real(v(1,:)); % force DC element real 0069 if rem(n,2) % odd output length 0070 x=real(ifft([v;conj(v(m:-1:2,:))],[],1)); % do it the long way 0071 else % even output length 0072 v(m,:)=real(v(m,:)); % force nyquist element real 0073 w=ones(1,k); 0074 % t=[cumprod([-0.5i; exp(2i*pi/n)*ones(m-2,1)]); 0.5i]; 0075 t=-0.5i* exp((2i*pi/n)*(0:m-1)).'; 0076 z=(t(:,w)+0.5).*(conj(flipud(v))-v)+v; 0077 z(m,:)=[]; 0078 zz=ifft(z,[],1); 0079 x=zeros(n,k); 0080 x(1:2:n,:)=real(zz); 0081 x(2:2:n,:)=imag(zz); 0082 end 0083 s(d)=n; % change output dimension 0084 if d==1 0085 x=reshape(x,s); 0086 else 0087 x=permute(reshape(x,s([d:ns 1:d-1])),[ns+2-d:ns 1:ns+1-d]); 0088 end 0089 end