## [636453]: fourier / dstiv.m  Maximize  Restore  History

### 83 lines (67 with data), 2.1 kB

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82``` ```function c=dstiv(f,L,dim) %DSTIV Discrete Sine Transform type IV % Usage: c=dstiv(f); % c=dstiv(f,L); % c=dstiv(f,[],dim); % c=dstiv(f,L,dim); % % DSTIV(f) computes the discrete sine transform of type IV of the % input signal f. If f is a matrix, then the transformation is applied to % each column. For N-D arrays, the transformation is applied to the first % dimension. % % DSTIV(f,L) zero-pads or truncates f to length L before doing the % transformation. % % DSTIV(f,[],dim) applies the transformation along dimension dim. % DSTIV(f,L,dim) does the same, but pads or truncates to length L. % % The transform is real (output is real if input is real) and % it is orthonormal. It is its own inverse. % % Let f be a signal of length _L and let c=DSTIV(f). Then % %M L-1 %M c(n+1) = sqrt(2/L) * sum f(m+1)*sin(pi*n*(m+.5)/L) %M m=0 %F \[ %F c\left(n+1\right)=\sqrt{\frac{2}{L}}\sum_{m=0}^{L-1}f\left(m+1\right)\sin\left(\frac{\pi}{L}\left(n+\frac{1}{2}\right)\left(m+\frac{1}{2}\right)\right) %F \] % % See also: dstii, dstiii, dctii % %R rayi90 wi94 % AUTHOR: Peter Soendergaard % TESTING: TEST_PUREFREQ % REFERENCE: REF_DSTIV error(nargchk(1,3,nargin)); if nargin<3 dim=[]; end; if nargin<2 L=[]; end; [f,L,Ls,W,dim,permutedsize,order]=assert_sigreshape_pre(f,L,dim,'DSTIV'); if ~isempty(L) f=postpad(f,L); end; s1=zeros(2*L,W); c=zeros(L,W); m1=1/sqrt(2)*exp(-(0:L-1)*pi*i/(2*L)).'; m2=-1/sqrt(2)*exp((1:L)*pi*i/(2*L)).'; for w=1:W s1(:,w)=[m1.*f(:,w);flipud(m2).*f(L:-1:1,w)]; end; s1=i*exp(-pi*i/(4*L))*fft(s1)/sqrt(2*L); % This could be done by a repmat instead. for w=1:W c(:,w)=s1(1:L,w).*m1+s1(2*L:-1:L+1,w).*m2; end; if isreal(f) c=real(c); end; c=assert_sigreshape_post(c,dim,permutedsize,order); % This is a slow, but convenient way of expressing the algorithm. %R=1/sqrt(2)*[diag(exp(-(0:L-1)*pi*i/(2*L)));... % flipud(diag(-exp((1:L)*pi*i/(2*L))))]; %c=i*(exp(-pi*i/(4*L))*R.'*fft(R*f)/sqrt(2*L)); ```