I'm trying to understand the behavior of two different methods of decimating and frequency translating a signal. At this point, I'm only swapping the translation and decimation; I'm not doing any polyrate processing. The graphic below depicts the basic idea. I have a linear phase filter around cf. The signal is translated by (for now offset is 0)
shift = -cf + FsOut/4 + offset
and then complex filtered, decimated down to FsOut/4, and converted back to real. The input rate is either 2, 4, 8, or 16 times the output rate.
I thought I might as well speed things up a bit by re-ordering the filter and decimation with the translation. Basically it's the same thing, except I multiply the frequency translation rate by the decimation factor. As it turns out, when offset is zero, the two different methods result in bit exact replicas of the signal. Note I'm quantizing to 8 bits. I believe when the decimation ratio is 16 and offset is a multiple of FsOut/8, the two methods may also be the same or at least somehow very close.
More generally, when offset is not zero, the time domain signal is out of phase, but I get approximately the same power spectrum. Now as frequency translation is not a linear operator, it's not immediately obvious why I would get the same time domain signal in this special case. I thought maybe it happens to work out that there is a simple linear relationship between the two different methods, so I spent some time trying to find a phase offset as a function of the offset and decimation ratio that would make the time domain signals always identical, but I couldn't come up with anything. I also tried looking at the phase of the cross-spectrum, which should be a line if one signal is just a delayed copy of the other, but it only works for offset equal to zero (obviously since they are identical).
After thinking about it for a while, my best guess is that if you were to write this as a poly-phase structure and commute the translation operation with the filter and decimate operations, it would just happen to work out that in the special case of FsOut/4, frequency translation could be pushed all the way across the final sum (after taking the real part). But normally, is it correct that there is no linear relationship between these two different methods, or have I just got something wrong here? Also, practically, is there any reason why it would matter which way I do it? I'm only comparing the two methods to make sure I've done it correctly.
I can't give out my code. But I can try to provide more detail. Also, I just realized I forgot to take the real part in method 2, but it doesn't appear to change anything. In method 1, the LPF (h) is centered about FsOut/4. So we have something like the following in Matlab.
f0 = shift/fsIn*2*pi;
y = x .* exp(1j*(0:blen-1)*f0);
y = filter(h,1,y);
y = real(y(1:dec:end));
In method 2, the same filter is located at cf. I reorder the operations.
y = filter(h,1,x);
y = y(1:dec:end);
f0 = shift/fsIn*2*pi*dec; % note decimation multiple
y = real(y .* exp(1j*(0:blen/dec-1)*f0));
To rephrase the question, is there a linear operator that takes the output of method 1 and yields the output of method 2 for offset not zero? I believe it only happens to work for FsOut/4 as described above. But I'm not sure...
Maybe there is a problem with translating the filter to high-frequency?
You were correct. Now that I've had some time to look this over, I see it was as simple as that. In fact, there is no correction factor needed at all. I don't know why I wasn't seeing a simple linear phase relationship between the two methods. I may have done something wrong when taking the cross-spectral phase, or maybe it was related to something else, like an issue with the quantization. Lesson learned. I should have distilled my simulation down to the simplest possible code when I got stuck as you did in your post.
So the issue was I was translating the upper filter to fc instead of (fc-offset). As you can see in the plot on top, the black lines are the center of the two filters. The red lines are where the filtering occurs. The center of the upper filter needed to be shifted by offset so that the phase change would be the same.
It's easy to calculate the phase offset introduced. The bottom diagram shows the linear phase slope for a unit delay filter. Using the filter I provided in the comments, the slope would be 640π/(FsIn/2). Thus the correction factor is exp(-i*640π/(FsIn/2)*offset). Of course, it doesn't make sense to fix the problem this way since you will still get some error due to ripple in the pass band. If you filter in the exact same place relative to the center frequency of each filter, the only error should be numerical.
I also noted that you happen get very small error if shift is a multiple of FsIn/(group delay) (e.g. 1e6/640), even with the offset error. The reason is obvious now; plug into the exponential above and you get exp(-i*2π) == 1. The error is probably due to slight magnitude response differences in the pass band and/or possibly due to the stop band of the negative sine frequencies. One last thing I noted empirically; there is no phase change due to decimation and aliasing of the upper filtered signal -- obviously, or this wouldn't work without phase correction.