methods of computing fixed point atan2 on FPGA I am in need of computing atan2(x,y) on an FPGA with a continuous input/output stream of data. I managed to implement it using unrolled, pipelined CORDIC kernels, but to get the accuracy I need, I had to perform 32 iterations. This led to a pretty large amount of LUTs being devoted to this one task. I tried changing the flow to use partially unrolled CORDIC kernels, but then I needed a multiplied clock frequency to execute repeated loops while still maintaining a continuous input/output flow. With this, I could not meet timing.

So now I am reaching out for alternative ways of computing atan2(x,y).

I thought about using block-RAM lookup tables with interpolation, but since there are 2 variables I would need 2 dimensions of lookup tables, and this is very resource intensive in terms of block-RAM usage.

I then thought about using the fact that atan2(x,y) is related to atan(x/y) with quadrant adjustment. The problem with this is that x/y needs a true divide since y is not a constant, and divisions on FPGAs are very resource intensive.

Are there any more novel ways to implement atan2(x,y) on an FPGA that would result in lower LUT usage, but still provide good accuracy ?

You can use logarithms to get rid of the division. For $(x, y)$ in the first quadrant:

$$z = \log_2(y)-\log_2(x)\\
\text{atan2}(y, x) = \text{atan}(y/x) = \text{atan}(2^z)$$

Figure 1. Plot of $\text{atan}(2^z)$

You would need to approximate $\text{atan}(2^z)$ in range $-30 < z < 30$ to get your required accuracy of 1E-9. You can take advantage of the symmetry $\text{atan}(2^{-z}) = \frac{\pi}{2}-\text{atan}(2^z)$ or alternatively ensure that $(x, y)$ is in a known octant. To approximate $\log_2(a)$:

$$b = \text{floor}(\log_2(a))\\
c = \frac{a}{2^b}\\
\log_2(a) = b + \log_2(c)$$

$b$ can be calculated by finding the location of the most significant non-zero bit. $c$ can be calculated by a bit shift. You would need to approximate $\log_2(c)$ in range $1 \le c < 2$.

Figure 2. Plot of $\log_2(c)$

For your accuracy requirements, linear interpolation and uniform sampling, $2^{14} + 1 = 16385$ samples of $\log_2(c)$ and $30\times 2^{12} + 1 = 122881$ samples of $\text{atan}(2^z)$ for $0 < z < 30$ should suffice. The latter table is pretty large. With it, the error due to interpolation depends greatly on $z$:

Figure 3. $\text{atan}(2^z)$ approximation largest absolute error for different ranges of $z$ (horizontal axis) for different numbers of samples (32 to 8192) per unit interval of $z$. The largest absolute error for $0 \le z < 1$ (omitted) is slightly less than for $\text{floor}(\log_2(z)) = 0$.

The $\text{atan}(2^z)$ table can be split into multiple subtables that correspond to $0 \le z < 1$ and different $\text{floor}(\log_2(z))$ with $z \ge 1$, which is easy to calculate. The table lengths can be chosen as guided by Fig. 3. The within-subtable index can be calculated by a simple bit string manipulation. For your accuracy requirements the $\text{atan}(2^z)$ subtables will have a total of 29217 samples if you extend the range of $z$ to $0 \le z < 32$ for simplicity.

For later reference, here is the clunky Python script that I used to calculate the approximation errors:

from numpy import *
from math import *
N = 10
M = 20
x = array(range(N + 1))/double(N) + 1
y = empty(N + 1, double)
for i in range(N + 1):
    y[i] = log(x[i], 2)

maxErr = 0
for i in range(N):
    for j in range(M):
        a = y[i] + (y[i + 1] - y[i])*j/M
        if N*M < 1000:
            print str((i*M + j)/double(N*M) + 1) + ' ' + str(a)
        b = log((i*M + j)/double(N*M) + 1, 2)
        err = abs(a - b)
        if err > maxErr:
            maxErr = err

print maxErr

y2 = empty(N + 1, double)
for i in range(1, N):
    y2[i] = -1.0/16.0*y[i-1] + 9.0/8.0*y[i] - 1.0/16.0*y[i+1]

y2[0] = -1.0/16.0*log(-1.0/N + 1, 2) + 9.0/8.0*y[0] - 1.0/16.0*y[1]
y2[N] = -1.0/16.0*y[N-1] + 9.0/8.0*y[N] - 1.0/16.0*log((N+1.0)/N + 1, 2)

maxErr = 0
for i in range(N):
    for j in range(M):
        a = y2[i] + (y2[i + 1] - y2[i])*j/M
        b = log((i*M + j)/double(N*M) + 1, 2)
        if N*M < 1000:
            print a
        err = abs(a - b)
        if err > maxErr:
            maxErr = err

print maxErr

y2[0] = 15.0/16.0*y[0] + 1.0/8.0*y[1] - 1.0/16.0*y[2]
y2[N] = -1.0/16.0*y[N - 2] + 1.0/8.0*y[N - 1] + 15.0/16.0*y[N]

maxErr = 0
for i in range(N):
    for j in range(M):
        a = y2[i] + (y2[i + 1] - y2[i])*j/M
        b = log((i*M + j)/double(N*M) + 1, 2)
        if N*M < 1000:
            print str(a) + ' ' + str(b)
        err = abs(a - b)
        if err > maxErr:
            maxErr = err

print maxErr

P = 32
NN = 13
M = 8
for k in range(NN):
    N = 2**k
    x = array(range(N*P + 1))/double(N)
    y = empty((N*P + 1, NN), double)
    maxErr = zeros(P)
    for i in range(N*P + 1):
        y[i] = atan(2**x[i])

    for i in range(N*P):
        for j in range(M):
            a = y[i] + (y[i + 1] - y[i])*j/M
            b = atan(2**((i*M + j)/double(N*M)))
            err = abs(a - b)
            if (i*M + j > 0 and err > maxErr[int(i/N)]):
                maxErr[int(i/N)] = err

    print N
    for i in range(P):
        print str(i) + " " + str(maxErr[i])

The local maximum error from approximating a function $f(x)$ by linearly interpolating $\hat{f}(x)$ from samples of $f(x)$, taken by uniform sampling with sampling interval $\Delta x$, can be approximated analytically by:

$$\widehat{f}(x) - f(x) \approx (\Delta x)^2\lim_{\Delta x\rightarrow 0}\frac{\frac{f(x) + f(x + \Delta x)}{2} - f(x + \frac{\Delta x}{2})}{(\Delta x)^2} = \frac{(\Delta x)^2 f''(x)}{8},$$

where $f''(x)$ is the second derivative of $f(x)$ and $x$ is at a local maximum of the absolute error. With the above we get the approximations:

$$\widehat{\text{atan}}(2^z) - \text{atan}(2^z) \approx \frac{(\Delta z)^2 2^z(1 - 4^z)\ln(2)^2}{8(4^z + 1)^2},\\
\widehat{\log_2}(a) - \log_2(a) \approx \frac{-(\Delta a)^2}{8 a^2\ln(2)}.$$

Because the functions are concave and the samples match the function, the error is always to one direction. The local maximum absolute error could be halved if the sign of the error was made to alternate back and forth once every sampling interval. With linear interpolation, close to optimal results can be achieved by prefiltering each table by:

$$y[k] = \cases{\begin{array}{rrrrrl}&&b_0x[k]&\negthickspace\negthickspace\negthickspace+ b_1x[k+1]&\negthickspace\negthickspace\negthickspace+ b_2x[k+2]&\text{if } k = 0,\\
&c_1x[k-1]&\negthickspace\negthickspace\negthickspace+ c_0x[k]&\negthickspace\negthickspace\negthickspace+ c_1x[k+1]&&\text{if }0 < k < N,\\
b_2x[k-2]&\negthickspace\negthickspace\negthickspace+ b_1x[k-1]&\negthickspace\negthickspace\negthickspace+ b_0x[k]&&&\text{if } k = N,
\end{array}}$$

where $x$ and $y$ are the original and the filtered table both spanning $0 \le k \le N$ and the weights are $c_0 = \frac{9}{8}, c_1 = -\frac{1}{16}, b_0 = \frac{15}{16}, b_1 = \frac{1}{8}, b_2 = -\frac{1}{16}$. The end conditioning (first and last row in the above equation) reduces error at the ends of the table compared to using samples of the function outside of the table, because the first and the last sample need not be adjusted to reduce the error from interpolation between it and a sample just outside the table. Subtables with different sampling intervals should be prefiltered separately. The values of the weights $c_0, c_1$ were found by minimizing sequentially for increasing exponent $N$ the maximum absolute value  of the approximate error:

$$(\Delta x)^N\lim_{\Delta x\rightarrow 0}\frac{\left(c_1f(x - \Delta x) + c_0f(x) + c_1f(x + \Delta x)\right)(1-a) + \left(c_1f(x) + c_0f(x + \Delta x) + c_1f(x + 2 \Delta x)\right)a - f(x + a\Delta x)}{(\Delta x)^N}
=\left\{\begin{array}{ll}(c_0 + 2c_1 - 1)f(x) &\text{if } N = 0, \bigg| c_1 = \frac{1 - c_0}{2}\\
0&\text{if }N = 1,\\
\frac{1+a-a^2-c_0}{2}(\Delta x)^2 f''(x)&\text{if }N=2, \bigg|c_0 = \frac{9}{8}\end{array}\right.$$

for inter-sample interpolation positions $0 \le a < 1$, with a concave or convex function $f(x)$ (for example $f(x) = e^x$). With those weights solved, the values of the end conditioning weights $b_0, b_1, b_2$ were found by minimizing similarly the maximum absolute value of:

$$(\Delta x)^N\lim_{\Delta x\rightarrow 0}\frac{\left(b_0f(x) + b_1f(x + \Delta x) + b_2f(x + 2 \Delta x)\right)(1-a) + \left(c_1f(x) + c_0f(x + \Delta x) + c_1f(x + 2 \Delta x)\right)a - f(x + a\Delta x)}{(\Delta x)^N}
=\left\{\begin{array}{ll}\left(b_0 + b_1 + b_2 - 1 + a(1 - b_0 - b_1 - b_2)\right)f(x) &\text{if } N = 0, \bigg| b_2 = 1 - b_0 - b_1\\
(a-1)(2b_0+b_1-2)\Delta x f'(x)&\text{if }N = 1,\bigg|b_1=2-2b_0\\
\left(-\frac{1}{2}a^2 + \left(\frac{23}{16} - b_0\right)a + b_0 - 1\right)(\Delta x)^2f''(x)&\text{if }N=2, \bigg|b_0 = \frac{15}{16}\end{array}\right.$$

for $0 \le a < 1$. Use of the prefilter about halves the approximation error and is easier to do than full optimization of the tables.

Figure 4. Approximation error of $\log_2(a)$ from 11 samples, with and without prefilter and with and without end conditioning. Without end conditioning the prefilter has access to values of the function just outside of the table.

This article probably presents a very similar algorithm: R. Gutierrez, V. Torres, and J. Valls, “FPGA-implementation of atan(Y/X) based on logarithmic transformation and LUT-based techniques,” Journal of Systems Architecture, vol. 56, 2010. The abstract says their implementation beats previous CORDIC-based algorithms in speed and LUT-based algorithms in footprint size.


Is it possible to calculate/estimate the storage space a certain duration of .mp4 video files will take up? Is it possible to calculate/estimate the storage space a certain duration of .mp4 video files will take up? I have the resolution (eg. 4000 X 3000), the format (H.264 or H.265), the audio format (AAC) and the sample rate (41 000 Hz).

No.

The MPEG4 video codec allows for a large ranges (easily a dynamic range of 50 for "useful" quality) of compression ratios for the same raw video material. The same applies to audio compressors.

And of course, there's simply video material that is easy to compress, so using the same type of quantization, motion prediction and frame structure, some videos will compress way better than other.

In fact, in many cases, you ask your encoder to aim for a target average bitrate even, so that you don't have to estimate anything, but set the amount of bits per second. If you have a technical necessity to know the bitrate beforehand, that would be the way to go!


Wrong output for phase demodulation I tried to plot phase demodulated signal and original message signal but there is a significant difference between magnitudes of the message and demodulated signal. I have attached the code and the graph. Please help in pointing out the error.

Code

clear all; close all;
fm = 50;
fc = 1000;
fs = 4000;

t = 0:1/fs:0.1;
kp1 = pi/2;

m_t = 5*cos(2*pi*fm*t);

x_pm = cos(2*pi*fc*t + kp1*m_t);

x = pmdemod(x_pm,fc,fs,kp1);

figure(1);
plot(t,x_pm);
title('phase modulated signal');

figure(2);
plot(t,x,t,m_t);
title('message and demodulated signal');
legend('demodulated signal','message signal');

plot

You tell the phase demodulator that the signal has a maximum phase deviation of $\pi/2$, but then you actually use a signal with a phase deviation of $5\pi/2$ (due to the amplitude of the message signal). This doesn't make sense.

Note that even the most ideal phase demodulator cannot distinguish between $\cos(2\pi f_ct+\phi(t))$ and $\cos(2\pi f_ct + \phi(t)+2k\pi)$, $k\in\mathbb{Z}$.


What is the secrecy rate? I read several papers on physical layer security that explain several methods to maximize the secrecy rate. However, I didn't find a very clear definition about the secrecy rate. What is it? How is it computed?

Here is a pretty clear paper explaining secrecy rate and uses over fading channels.

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.105.2844&rep=rep1&type=pdf

See equation 3 for the formal definition. One way to think about it is like a tug-of-war between the main channel and wiretapping channel. If the main channel has a higher capacity then the secrecy rate is positive, otherwise it is negative.


How do I reconstruct text from an image using only morphological operations? I want to reconstruct the text from the following image as best as possible. The tricky part is that I want to do it only using morphorogical operations on the image

I tried using erosion, dilation, opening and closing but the result is not very good.

Is this even possible?

If you're willing to add/subtract etc. morphologically transformed images, you can count how many signal pixels are in the vicinity of each pixel, and threshold based on that number.

img = imread('http://i.stack.imgur.com/wicpc.png');

n = false(3);n(4) = 1;
s = false(3);s(6) = 1;
w = false(3);w(2) = 1;
e = false(3);e(8) = 1;

%# note that you could convolve with a cross instead
fourNeighbourCount = imerode(img,n) + imerode(img,s) + imerode(img,w) + imerode(img,3) + img;

%# require at least two neighbours
img = fourNeighbourCount > 1;

If you then convolve with a 3x3 mask that has a hole in the middle, you can get something like this:


Intuition behind cepstrum analysis terminology In multiple books I found a parallel between spectrum-analysis terminology and cepstrum-analysis terminology that says

 low pass filter $\rightarrow$ short pass lifter

I didn't read explicitly what lifter does, but assuming a short pass lifter passes short rahmonics just as low pass filter passes low harmonics, why the parallel isn't rather

 low pass filter $\rightarrow$ long pass lifter

since low harmonics correspond to long rahmonics?

in filtering we try to separate two signals from each other. if two signals are added together this is the normal filtering. but these two signals could be multiplied or convoluted with each other. considering there are huge scientific background on this kind of filtering ( for additive signals) we try to transform other problems to this kind.

if signals multiplied with each other, we take the logarithm of their multiplication to convert it two additive signal. then we use filtering methods for additive signals.

if they were convoluted with each other, we know their spectrum is multiplied, so by taking the logarithm of their spectrum, we reach two additive spectrum in frequency domain. now we treat these two additive spectrum in frequency domain as two additive signals and apply the filtering methods for additive signals. considering we work on frequency domain signals if we use normal filtering terminology there would be some confusion e.g. if we say we use a low-pass filter someone might think we apply a LPF on our original time domain signal instead of their frequency domain counter part, so in order to avoid these confusions some new terminology is used.

now if we take Fourier transform of these frequency domain signals we go to time domain, so a low-pass filter (for frequency domain signals) becomes a short-pass filter which passes the signal at times before a cut-off time and high-pass filter becomes long-pass filter which passes the signals after a cut-off time.


Speech signal formant spectrogram. How can we determine F1? The formant spectrogram in this page points F1 and F2 as below;

In bottom-right F1 movement, I thought F1 could be below blue one. But how can we determine F1 is upper one not below one?

Is it because we know F1 and F2 have to be getting closer from [ε] to [ae]? Or, can we determine F1 only with spectrogram?

You should be able to see formants very clearly as peaks on the spectrogram slices. The 1st, 2nd ... are the several peaks. In the picture above (in the blue area) the first two peaks seem to have very similar power and so this can introduce errors in the estimation algorithm.

You can use Praat and "Show formants" option to track formants and observe their behaviour over time. You can try to make continuous transitions from vowel to vowel and then try the same vowels with a silence in the middle and observe the differences.

If you are developing an algorithm you can use a Kalman filter for tracking formants. (maybe this helps: http://link.springer.com/content/pdf/10.1007/978-3-540-85980-2.pdf#page=257)


What is a sliding window? I tried to find the answer by my own but unfortunately I couldn't.

Let's say I have the discret values of a signal which is 10s long.

What is the FFT which is performed on a sliding window with time?

What is the difference to the FFT, where I get the frequency spectrum of the signal?

If you have the samples of a signal of length $L$ (corresponding to 10s analog duration for example) and want to analyse its spectrum, then you can peform an $L$ point DFT/FFT on the whole data at once.

This would provide you the maximum spectral resolution but minimum temporal resolution; i.e., time locality of the events are mainly lost (or better become obscured in the phase spectrum).

If you wish to sacrifice some spectral resolution but increase the temporal instead, then you can analyse the signal in small blocks (overlapping chunks) and move the block window in time to observe the change in the local spectral behaviour. This motion of the analysis window is referred to as sliding action. And the analysis is essentially a sliding window analysis, aka short-time Fourier analysis, windowed Fourier analysis etc.

Note that for various reasons it's best to use some weighting within windowing before the DFTs are applied.


What segmentation methods can be used for simple images? As related to the project I am currently involved in, I have come to realize that one thing that might help is for me to 'segment' my image. To that end, I have shown a very simple example below.

How do segmentation methods work in general, and what class of them would help getting me from the first image:

...to the second:

Please ignore the scaling that might seem to have occured here, but you get the general idea.

One of the simplest method you can use is

"Region growing" scheme. Here the algorithm starts with a seed points, and keeps registering the pixel surrounded by it as long as the region as a whole has same property as before and after the inclusion of pixel. This is applied with reasonably good set of seed points recursively.

This property can be color variations, statistical properties such as mean and variance of intensity of the regions, or it can be texture. And the criteria for segmentation is homogeneity of the region.

The same principle but alternative implementation is region split and merge approach. Here a image is first divided into fixed parts - say 4 parts and corresponding scores of above properties is computed. Now the each sub region is split further, and we can check if any region is significantly different from others or same. If homongenety of sub parts is same essentially it is merged back. But if the element is realized to have different properties it now fall in different regions.

Read this for more details: http://vplab.iitm.ac.in/courses/CV_DIP/PDF/lect-Segmen.pdf


How do you alter the phase relationship between I and Q in a quadrature signal? I am collection a quadrature signal from a radio.  The I and Q signals go through a sound card at a sample rate of 192000, and then a FFT is performed on the signals.
Due to various reasons, the signals are shifted slightly from the 90 degree phase difference and need to be corrected at some point. How can one go about doing this?
Thanks for your help.
Tom

As you already know, the quadrature error, also called quadrature skew, describes how far off the actual angle between I and Q is from the ideal 90 degrees.
It is one component of modulation error ratio.

There seem to be lots of people talking about how to measure quadrature skew,
but few people talking about how to compensate.
The two major approaches to compensation seem to be:

During initial factory tuning, some human watches an o'scope measuring the skew while manually tweaking the analog parts, or
During normal operation, the receiver internally estimates the quadrature skew and compensates digitally.

Often hardware designers deliberately build "off-tuned recievers".
With such hardware, when the transmitter repeatedly transmits "the same" constellation point, the raw I/Q point sampled by the reciever -- rather than staying close to a single fixed location in the I/Q diagram as with a perfectly tuned direct-conversion receiver -- instead rotates in a circle.

Such a receiver somehow estimates (typically in a FPGA or in software) how far the constellation has been rotated by that analog hardware, and corrects for it by rotating the I/Q sample in the opposite direction.
It is possible to simultaneously estimate the phase skew between I and Q and correct for that at the same time.
(In effect,
if you know the skew is S and you measure raw x,q coordinate,
you can deskew and find the true x,y in perfectly orthogonal coordinates
with y = qc - xt,
where c is some approximation of 1/cos(S) and t is some approximation of tan(S); often c = 1 and t=S (in radians) is adequate).

Further reading:

"Blind Equalization of Receiver In-Phase/Quadrature Skew in the Presence of Nyquist Filtering"
a patent search for "quadrature skew" compensation turns up a few.
"Effects of physical layer impairments on OFDM systems"
"Compensation for In-Phase/Quadrature Imbalance in Coherent-Receiver Front End for Optical Quadrature Amplitude Modulation"
"Correcting Imperfections in IQ Modulators to Improve RF Signal Fidelity"
"Estimation and Compensation of Impairments caused by Quadrature Modulators and Demodulators"


What is the smallest sampling points that can recover the signal? Propose that we have a signal with a model of
$$
f_{4at0} = v_0(1 - \epsilon + \epsilon \exp(-c_l(t_{t2d}-t_0)))
$$

$f_{4at0}$ is dependent variable while $t_{t2d}$is independent variable,others are parameters.

It is a representation of the relationship of biological material concentration-time that is inside the blood serum.Since this material cannot be assayed real time in the human body, the blood samples have to be collected from patients.It's really inconvenient to obtain blood samples from human.So, my question is what's the lower limit of sampling points? And how to schedule the sampling points?

ps: here is the matlab code that I used .I want to use the "Nyquist Sampling Principle",but I don't know how.

cl=6.2;

epsilon=0.95;

t0=0.3625;

delta=0.1;

v0=1000000;

f4at0=v0*(1-epsilon+epsilon*exp(-cl*(tt2d-t0)));

syms tt2d;

f4=v0*(1-epsilon+epsilon*exp(-cl*(tt2d-t0)));

fourier(f4,tt2d)

Output:

  ans =

  100000*pi*dirac(-tt2d) + 950000*transform::fourier(exp((31*tt2d)/5 + 899/400), -tt2d, -tt2d)

I don't know how to understand transform::fourier(exp((31*tt2d)/5 + 899/400), -tt2d, -tt2d).

I don't know how to calculate T since there is also a Dirac function.
I don't know if there are other principles apart from the "Nyquist Sampling Principle" that can be used here.

This depends if you want to do sampling or just model identification. Model identification means "figuring out the numbers in your model from actually measured data". If the model is good and the data is not noisy you only need a few points. The minimum is equal to the number of unknown model parameters. Let's say you don't know v0, cl, epsilon and t0. In this case you need at least four points to identify those. Your model isn't linear so there is some math involved, however it's not overly complicated and very workable.
If the measurements are noisy it would be good to take more points and do a least square error fit. The number of points depends on the amount of noise and precision required for the results.

You can also sample the sequence but you need to consider sampling theorem. Your time sequence represents a first order lowpass filter. The cutoff frequency is determined only through your model parameter "cl" and is simply cl divided by 2 pi. At this frequency the spectrum is at -3 dB, at ten times this frequency it's at -20dB and at a 100 times this frequency it's -40 dB etc. Since it never gets all the way to zero you have to live with some amount of aliasing. Again it boils down to how much error you can/want to tolerate.


How to create a sine wave generator that can smoothly transition between frequencies I am able to write a basic sine wave generator for audio, but I want it to be able to smoothly transition from one frequency to another. If I just stop generating one frequency and immediately switch to another there will be a discontinuity in the signal and a "click" will be heard.

My question is, what is a good algorithm to generate a wave that starts at, say 250Hz, and then transitions to 300Hz, without introducing any clicks. If the algorithm includes an optional glide/portamento time, then so much the better.

I can think of a few possible approaches such as oversampling followed by a low pass filter, or maybe using a wavetable, but I am sure this is a common enough problem that there is a standard way of tackling it.

One approach that I have used in the past is to maintain a phase accumulator which is used as an index into a waveform lookup table. A phase delta value is added to the accumulator at each sample interval:

phase_index += phase_delta

To change frequency you change the phase delta that is added to the phase accumulator at each sample, e.g.

phase_delta = N * f / Fs

where:

phase_delta is the number of LUT samples to increment
freq is the desired output frequency
Fs is the sample rate

This guarantees that the output waveform is continuous even if you change phase_delta dynamically, e.g. for frequency changes, FM, etc.

For smoother changes in frequency (portamento) you can ramp the phase_delta value between its old value and new value over a suitable number of samples intervals rather than just changing it instantaneously.

Note that phase_index and phase_delta both have an integer and a fractional component, i.e. they need to be floating point or fixed point. The integer part of phase_index (modulo table size) is used as an index into the waveform LUT, and the fractional part may optionally be used for interpolation between adjacent LUT values for higher quality output and/or smaller LUT size.


Meaning of "Base" in Walsh Hadamard Tranform in Image Processing? Currently, I am studying a paper by K. He and J. Sun, “Computing nearest-neighbor ﬁelds via propagation-assisted KD-trees,” in Proc. IEEE Conf. Comput. Vis. Pattern Recog., 2012, pp. 111–118.

In the paper there is a sentence which i can quote as:

  ...  We use the ﬁrst 16 WHT bases for the Y channel2 and 4 for each
  chrominance channel (Cb/Cr) throughout this paper. ...

Although I have made a search on the internet, I couldn't get the meaning of base.

In MATLAB 2014b There is a function fwht(x,n) where n is the order of the transform. Can this order be the base that i am trying to find the meaning of?

I believe that here bases stands for "basis vectors". From the cited paper by  Y. Hel-Or and H. Hel-Or. Real time pattern matching using projection kernels, 2003, ICCV,  you can see the set of $8 \times 8$ vectors from  a separable Walsh-Hadamard basis:

[20] Y. Hel-Or and H. Hel-Or. Real time pattern matching using
projection kernels. In ICCV, pages 1430–1445, 2003.

It consists in a tensor product of two $8$ 1D Walsh-Hadamard transforms (one horizontal, one vertical).

The sixteen 2D vectors  in the left-top $4\times 4$ square are the one with the lowest sequencies (sequency is the  frequency analog for Walsh-Paley-Hadamard transforms), in terms of "maximum sequency":


what are and why are sine and cosine modulated integrals used? I have found the definition of the following formulas in a paper regarding active vibration control, where they are called sine and cosine modulated integrals.

$y$ is measurement signal with a strong periodic component of frequency $N\Omega$

$$y^{(i)}_{Nc}=\frac{2}{T}\int^{T}_{0}y^{(i)}(ϕ)\cos(Nϕ)dϕ$$

$$y^{(i)}_{Ns}=\frac{2}{T}\int^{T}_{0}y^{(i)}(ϕ)\sin(Nϕ)dϕ$$

where $\phi=\Omega t$.

From these the vector $y_N^{(i)}$ is defined as
$y_N^{(i)} = \begin{bmatrix} y_{Nc}^{(1)}\\ y_{Ns}^{(1)}\\\vdots\end{bmatrix}$

The same is done for the control input(s) $u$.
Then a quadratic cost function to be minimised at each step is defined using these newly introduced signals in this way:

$J(k) = y^T_NQy_N+u^T_NRu_N $ where $Q$ and $R$ are just two weighing matrices.

Are they introduced to consider the sine and cosine components at the considered frequency of the signal?

Can someone explain their meaning and why are they introduced?
Here the link for the paper

Another question: suppose I have the value $y_N$: how can I invert the relationship to get $y^{(i)}$

Here the link for the paper.

As someone made me notice in the comments, the integrals are just the coefficients of the Fourier Series for that frequency component (when the series is expressed with the cosine and sine parts separately).

The integral is performed on a period much bigger then the period corresponding to the frequency of the disturbance so that a better averaged is obtained.

To get back the signal in the time domain it is sufficient to multiply the two components by the cosine and sine respectively. Of course in this way the signal has only one frequency component.


Diode demodulation for AM What is the capacitor and resistance at output of LPF meant for in the document it was mentioned that it is DC blocker I don't get it can't message contain DC value.

Because of the diode at the input, you get a DC offset at the output. Notice that a single diode is already a primitive rectifier, because it blocks the negative half wave.

But you usually will want a purely AC output, that's why you put a HPF with very low cutoff frequency after the LPF. It will block the DC and yield a pure AC output signal. You can see it either as part of the demodulator or as part of the following circuit, usually some kind of amplifier or buffer, but it is always required.


How do you measure "detail" of a signal? I have an image and I would like to measure the amount of detail in it. Another way to look at it is to measure how blurry an image is. One way is to analyse the high frequency components in the Fourier transform of the image.

Are there any other/better methods?

What you are referring to is typically known as "Image Sharpness". A quick scan, as well as some prior knowledge, comes to the following:

Fourier analysis- Using this has 2 key disadvantages. First of all, noise would tend to show up no matter what, and thus higher frequency components would tend to show up. Secondly, sharpness tends to be a local phenomena, and thus might not show up if you do a transform of the entire image.
Eigenvalue analysis- I haven't actually read this paper, but it proposes using eigenvalue analysis to determine the sharpness of an image.
Edge detection algorithms depend on a certain amount of sharpness. One could use different values for edge detection parameters to determine the amount of sharpness.
Kurtosis Measurement of wavelet coefficients- Again, I haven't read the entire paper, but this seems to suggest calculating wavelet coefficients, performing an FFT of the entire set of coefficients, and measuring the kurtosis. This should be relatively immune to noise.

I'm sure there are many more. This is a very active field of study currently. If none of these methods suits you, then continue to search through academic papers, and see if you can find a better method.


Obtaining the Impulse response from the frequency of a discrete system I have the following discrete system's transfer function (Z domain):

$$h[z] = {z \over z - \frac{1}{2} }$$

I need to obtain the following:

The frequency response.
The impulse response.
The fourier spectre.

The first one is pretty easy we just substitute $z$ by $ e^{j \omega}$ giving us

$$h[\Omega] = { e^{j \omega} \over  e^{j \omega} - \frac{1}{2} }$$

In the second one is where I have some questions, as a rule the IR is just the anti-transform of the transfer function $h[\omega]$, the original transfer function $h[z]$ maps quite easily to $h[n] = ({1 \over 2})^n u[n]$, however I read or heard, not sure where though that the impulse response obtained from the frequecy response must be continuous, under that assumption we would replace $e^{j \Omega}$ with $\$$ giving us $h(\$) = {\$ \over $ - {1 \over 2}}$, the we would Laplace anti-trasnform to obtain the impulse response, so which approach would be the correct one? maybe I'm confusing concepts.

Now for the sake of completenes for the fourier espectre we just obtain the magnitude ($|h[\omega]|$) and phase of $\theta_{h[\omega]}$ and evaluate them over a set of $\omega$ values, plotting the results.

Thanks in advance.

As the frequency response of the system is asked, the RoC woulbe be $|z|>\frac{1}{2}$. Thus, the impulse response is right sided and is given by $$h[n] = \left(\frac{1}{2}\right)^n u[n]$$ We obtain the frequency response by replacing $z = e^{j\omega}$ in $H(z)$. Thus, $$H(e^{j\omega}) = \dfrac{1}{1-\frac{1}{2}e^{-j\omega}}$$


Applying a lowpass filter to a noisy square signal leads to a shift of the signal My knowledge of analog and digital signal analysis is behind me.. thus I do not have enough confidence in my procedure to analyze my signal, and I do not know which direction to pursue to improve. I'm measuring a noisy square signal at a sampling frequency fs = 1e6 Hz. I can increase this frequency up to 1e9 Hz if needed.
At a frequency fs = 1e6 Hz, I get the following signal:

Above, you can see 2 examples with 2 square waveforms at 2 different amplitudes. As you can see, the signal is prone to overshoot, in the example below, before the square starts; and is also noisy. Thus, I wanted to clean the signal by applying a low-pass filter.
I tried the 2 following: a Butterworth and a bessel filter using the scipy library in Python.
def butter_lowpass_filter(data, cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    y = lfilter(b, a, data)
    return y

def bessel_low_pass_filter(data, cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = bessel(order, normal_cutoff, btype='low', analog=False)
    y = lfilter(b, a, data)
    return y

I then call both filters on my data with the settings cutoff = 100000 # Hz, fs = 1e6 # Hz, order = 8. Why those settings? That is my first problem, I don't really know how to choose the cutoff frequency or the order. I simply tried different values, and those seem to be nice.
The result looks like this (left, low amplitude pulse; right, higher amplitude pulse):

The bessel filter seems to work better than the Butterworth, however both have the same problem. The signal is shifted to the right. If I zoom , I observe the clear shift:

The opposite shift seems to be present at the end of the square. Is it possible to correct this shift?
My end goal is to measure both the amplitude and the width of the square. Do you think my approach is correct or would you recommend me another approach?
Thanks for the guidance,
Mathieu

The only possible way that any real-world filter can work is by using past values of the signal you're filtering to get the latest filtered value.  So any real-world filter will always have some delay.
So if by "correct this shift" you mean "make it go away" -- not really.
However, if you want to look at the signal and its filtered version after the fact, then yes -- calculate the filter's group delay (or just eyeball it in your data), and in your graphs shift the input over by that amount.  This is easier to do if you use a symmetric FIR filter, because it'll have a constant group delay and the group delay is just half the filter length.
If your goal is to measure the duration and magnitude, then as long as you can have some delay in the measurement then just take your filtered version and feed it into whatever duration-and-magnitude measuring algorithm you've cooked up.  Since the leading and trailing edges of the square wave are offset by the same amount, the measured duration should be just fine.


How to choose between taking the real part or the absolute value of an inverse discrete Fourier transform? I kind of expected this to of been asked before, so if this is duplicate please point me in the right direction, however I couldn't find anything.

I am currently experimenting with low pass filtering. I take an image and apply a 2-dimensional fast Fourier transform to the image, next I apply a shift such that the zero frequency component is moved to the center (Using Matlab this is fftshift). Then I zero Fourier co-efficients which are not within some distance $r$ of the center pixel. I apply the inverse of the shift and then an inverse 2-dimensional fast Fourier transform.

Now I end up with a matrix of complex values, but to display an image I want a matrix of real values. With the matrix obtained after applying the IDFT I mapped the pixel values to their real part, imaginary part, magnitude and phase so that I could view the matrix as an image. I found that the real part and the magnitude looked like the original image and taking the imaginary part and phase rendered the image unrecognizeable compared to the original image. I also found that their wasn't a recognizable difference in quality between taking the real value of the transformed image or the absolute value.

My question is, is there advantages or disadvantages between taking the real part of absolute value after filtering (low pass, or in general)?

This is my first question on DSP.SE, so please leave comments on how to improve my question if it's not up to standard

If your processed 2D FFT is conjugate symmetric, then any imaginary components in the complex IFFT result are just numerical noise or rounding errors.  Thus using them to compute an absolute value is mostly a waste of CPU cycles.


How does a wavelet help in image compression? I understand how the discrete cosine transform is used in image compression in standards like JPEG. However, the concept of wavelets is a mystery to me. I do know that wavelets are functions that have a wave like functions with value >0 only for short intervals of time, then they vanish. That is why they are called wavelets. I do not know anymore

What is not clear to me is precisely how such functions help in image compression and achieve a higher compression rate than our other transforms like the cosine transform.

What is the best place to learn this?

A common wavelet based standard is JPEG 2000 and a common DCT based standard is JPEG.

JPEG 2000 uses wavelets, but a good portion of the better compression it achieves than JPEG is due to the fact that JPEG uses a much much simpler entropy coder (JPEG does context-dependent Huffman codes and run length coding, JPEG 2000 does arithmetic coding with some extra tricks). You can gain a lot by tweaking the entropy coding stage. The particular wavelet transform you choose will make a difference as well. And how you decide to compress the coefficients makes a huge difference as well (JPEG reads its dct block coefficients in a zig-zag pattern and uses run length coding to compress long runs of zeros in the higher freq components due to the nature of the quantizers it uses).

The primary advantage of wavelets as used in JPEG 2000 versus rounding the DCT coefficients in the manner of JPEG is that it reduces blocking artifacts.

If you're not familiar with wavelets, you can look at an introductory text on wavelets, like Mallat's Wavelet tour of Signal Processing or Vetterli's new book (freely available online). Then, you can read the JPEG 2000 standard. These notes are kinda nice at a high level as well (and it has good pointers to the references).

For a quick review of JPEG, you can look at Gonzalez & Woods' Digital Image Processing (second or third edition should be fine).


Designing observer with non-observable system I have this state-space system

\begin{align}
\dot{x}&=\begin{bmatrix}1 & 0\\3 & -2\end{bmatrix}x+\begin{bmatrix}10\\0 \end{bmatrix}u\\
y&=\begin{bmatrix}1 & 0\end{bmatrix}x
\end{align}

And I'm asked to design an observer such that its poles are at $-1+j$ and $-1-j$.

The problem I'm having is that the system described above is not observable, as its observability matrix's rank is 1.

Is there anything I can do in a case like this?

If the unobservable dynamics is stable, you should be able to separate the unobservable states and the observable states via a similarity transform. You then design an observer for the observable part.


Partial column approximation error Can someone help me understanf how to plot the next partial column in MATLAB:

while N=2001

|M|<700

ak =

N = 2001;
M = 600 ;
vectorOfK = (-M:M)';

vectorOfN = (0:N-1);
%%% For a given n
% n  = 314;

vectorOfAk = sin(199*pi*vectorOfK/2001)./(2001*sin(pi*vectorOfK/2001));
% Modify the NaN value (0/0 division) at k = 0;
vectorOfAk(M+1) = 199/2001;

matrixOfExpo = exp(1i*2*pi/N*(vectorOfK*vectorOfN));
matrixOfAkExpoProduct = repmat(vectorOfAk,1,N).*matrixOfExpo;

vectorOfATilde = sum(matrixOfAkExpoProduct)/N;

subplot(2,1,1)
plot(vectorOfN,real(vectorOfATilde));
xlabel('Real part (index)');xlabel('Amplitude (a. u.)');

subplot(2,1,2)
plot(vectorOfN,imag(vectorOfATilde));
xlabel('Imaginary part (index)');xlabel('Amplitude (a. u.)');


Why is the number of frequencies decomposed in scipy.signal.stft() equal to the hop size? I asked this over on StackOverflow first but figured you fine people might know more about signal processing (duh).

I'm trying to calculate the spectrogram of an audio file in python. There's a nice library function for doing this in scipy.signal.stft(), but it's not behaving as I expect it to.

Since a spectrogram is just a Short-time Fourier Transform, I expect the library function to chop my audio file up into segments ('time frames'), and then do a Fourier transform on each time frame.

This is what the function does as far as I can tell, but the size of the output matrix is behaving strangely. I'll explain:

The parameters you can specify for the function is nperseg and noverlap, which are the number of audio samples in each time frame and how much each time frame should overlap with the previous respectively. From this, we can caluclate the hop size, i.e. how many samples should be between each time frame: H = nperseg - noverlap.

This is where the strangeness comes in. The number of frequencies 'measured' in the Fourier transform of each time frame is exactly equal to the hop size, and it's always evenly spaced between 0 Hz and half the audio sampling frequency.

Say I have an audio file sampled at 22 kHz and I've set the function parameters so that the hop size is, say, 5 audio samples between each time frame. Then the output will show 5 different frequencies evenly spaced between 0-11 kHz, like this: [0kHz  2.75kHz  5.5kHz  8.75kHz  11kHz].

Why only these 5? What if I want to have a hop size of 5 but also see what the Fourier transform looks like for more frequencies? What if I want to have a hop size of 706 samples and 149 measured frequencies between 50 kHz and 65 kHz? Is this impossible?

The number of frequencies 'measured' in the Fourier transform of each time frame is exactly equal to the hop size,

The number of frequencies is not determined by the hop size, it's determined by the FFT size.  If you don't specify the FFT size, it uses nperseg as the FFT size.

  and it's always evenly spaced between 0 Hz

Yes, FFT is always evenly spaced, and always starts at 0 Hz.  It doesn't implement Zoom FFT (since it's meant primarily for reconstruction with ISTFT, and Zoom FFT would throw away the required information).

  and half the audio sampling frequency.

Yes, your input is real, so the negative frequencies are redundant, so they are not returned.   nperseg is 100, so nfft is 100, so it returns 100//2+1=51 frequencies (including 0 and Nyquist).

If you give it a complex input, then len(f) will equal nperseg because it returns both the positive and negative frequencies.

  What if I want to have a hop size of 706 samples and 149 measured frequencies between 50 kHz and 65 kHz?

It doesn't support Zoom FFT, so, assuming your sample rate is 65*2 = 130 kHz, you'd need to calculate how many samples you want across the entire bandwidth (= 1292?) and enter that as nfft, and specify noverlap = nperseg - hop, and then ignore the frequencies below 50 kHz that are returned.

Note that specifying nfft > nperseg returns redundant interpolated information that isn't necessary for reconstruction.


Checking convolution property II of z transform I have two sequences x and y of lengths, say 5 and 10.  I multiply them in time domain element by element.  I get another sequence.  Now this as per the convolution theorem should be equal to the inverse z tranform of the convolution of the two z transforms of the original sequences.  Now, I know how to get the ZTs of the original sequences.  How do I get the convolution of  these ZTs?

The term-by-term product of the sequences $(a_0, a_1, a_2)$ and
$(b_0, b_1, b_2)$ is $(c_0, c_1, c_2) = (a_0b_0, a_1b_1, a_2b_2)$. The corresponding $z$-transforms are
\begin{align}
A(z) &= a_0 + a_1z^{-1} + a_2z^{-2}\\
B(z) &= b_0 + b_1z^{-1} + b_2z^{-2}\\
C(z) &= a_0b_0 + a_1b_1z^{-1} + a_2b_2z^{-2}
\end{align}
According to the formula in Wikipedia pointed out to you by MattL, $C(z)$ can be obtained from $A(z)$ and $B(z)$ via a contour integral in
the complex plane.  The integrand of the contour integral (which
is with respect to the complex variable $v$) is
\begin{align}
v^{-1}A(v)B\left(\frac zv\right)
&= v^{-1}\left(a_0 + a_1v^{-1} + a_2v^{-2}\right)
\left(b_0 + b_1\left(\frac zv\right)^{-1}
+ b_2\left(\frac zv\right)^{-2}\right)\\
&= \frac {a_0b_0 + a_1b_1z^{-1} + a_2b_2z^{-2}}{v}
+ \sum_{m=0}^2\sum_{n=0, n\neq m}^2
\left(a_mb_nz^{-n}\right) v^{-m + n -1}\tag{1}
\end{align}
where in that last sum in $(1)$, the exponent of $v$ in each term is either $\geq 0$ or $\leq -2$. In other words, if we regard each of the $8$ terms in that sum as a separate function of $v$, then the function
either has no poles inside the unit circle ($-m+n-1 \geq 0$) or
has a pole of multiplicity
greater than $1$ at the origin ($-m+n-1 \leq -2$).
Now, write the contour integral of the right side of
$(1)$ as the sum of $9$ different contour integrals and
apply Cauchy's integral theorem to evaluate each of these $9$ contour integrals. The conclusion is that the integrals of all the terms
in that sum in $(1)$ integrate to $0$ while the contour integral
of $\displaystyle \frac {a_0b_0 + a_1b_1z^{-1} + a_2b_2z^{-2}}{v}$
(which has a single pole at the origin) gives the desired
$C(z) = a_0b_0 + a_1b_1z^{-1} + a_2b_2z^{-2}$.
Generalization of all this to sequences of arbitrary lengths is
straightforward.

If you don't know
enough complex-variable theory to figure out why it is that the contour
integral of the first term in $(1)$ is just the numerator and why the integral of the other terms is $0$, then you
will not be able to carry out the verification that you desire.


noise floor of attenuator NF of attenuator is same as the attenuation value as the noise floor remains same but signal level gets attenuated? But NF of amplifier depends on the device as it increases both signal and noise?

The Noise Figure of a passive attenuator is the attenuation value of the attenuator. This is because noise figure is by definition the difference in the SNR of the output and the SNR of the input, with the input terminated with the characteristic impedance of the system (typically 50 ohms).  Using this definition, the noise figure of an amplifier similarly is the additional increase in noise beyond the gain of the amplifier when the input of the amplifier is terminated.
The ever-present thermal noise floor is a white noise source with a power spectral density given by kT with k as Boltzmann’s Constant and T as temperature. At room temperature in power units of dBm (dB relative to 1 mW) this ends up being a power spectral density of -174 dBm/Hz. In the case of an attenuator, thermal noise is the only noise source given it is a passive device. If the attenuator is matched to 50 ohms, and with the input tied to a 50 ohm source resistor only and no other gain at the input, the thermal noise at the input will be -174 dBm/Hz (kTB) at room temperature.
The thermal noise at the output of the 50 ohm attenuator will also still be -174 dBm/Hz! The attenuator cannot somehow attenuate the thermal noise that is present in all parts of the system-- it certainly attenuates any noise that is applied to the input, but in terms of power density, the same level of thermal noise is present at the output as a local independent noise source.
Thus in terms of (SNR out) - (SNR in): Consider a noise-free signal at the input with it's noise floor then limited by thermal noise alone; this sets the input SNR. The signal will certainly be attenuated by the attenuation value, but as we explained above the noise can't get pushed down lower than its theoretical limit (the thermal noise floor). Since the noise is not affected as it is the same value at the output and input in this case, the SNR is therefore changed by the attenuator value, which by definition is the Noise Figure.
Cascaded Noise Figure is then considering what occurs when you place an amplifier in front of an attenuator. Now in this case amplified noise is applied to the input, not the thermal noise floor specifically. This amplified noise like any other waveform can of course be attenuated, down to the limit of the thermal noise floor. Consider these two cases and how the overall Noise Figure compares:
The formula for cascaded Noise Figure is as follows below for two components, expandable for any number of components. I'll give the equation first but then demonstrate with two cases where the explanation above becomes very clear.
$$NF_{cascade} = NF_{N-1} + \frac{NF_N-1}{G_{N-1}}$$
This states that the overall noise figure of N cascaded devices is the noise figure of the first component plus the noise figure of the subsequent components reduced by the cumulative gain of to each components input. The gain (G) and noise figure (NF) are given in power ratio (non-dB) quantities.
Consider this with these two components:
Amplifier with 20 dB gain, and 3 dB Noise Figure
10 dB Attenuator
Case 1: Amplifier followed by Attenuator
According to the formula:
$NF_1= 10^{3/10} = 2$
$NF_2 = 10^{10/10} = 10$
$G_1 = 10^{20/10} = 100$
$NF_{cascade} = 2 + \frac{10-1}{100} = 2.09$
The Noise Figure in dB is $10 \log_{10}(2.09) = 3.2 \text{ dB} $
Notice the attenuator had very little effect due to the amplifier in front of it! This is the motivation for Low Noise Amplifiers (LNA) at the front-end of the receiver with very low loss components in front of it (typically a necessary filter).
Case 2: Attenuator followed by Amplifier
$NF_1= 10^{10/10} = 10$
$NF_2 = 10^{3/10} = 2$
$G_1 = 10^{-10/10} = 0.1$
$NF_{cascade} = 10 + \frac{2-1}{0.1} = 20$
The Noise Figure in dB is $10 \log_{10}(20) = 13 \text{ dB} $
Intuition:
Consider a received signal that has a 30 dB SNR, with the noise limited by thermal noise alone. In case 1, we amplify both the signal and the noise by the amplifier with 20 dB gain. This amplifier has a 3 dB noise figure which means the signal goes up by 20 dB, but the noise goes up by 23 dB! (Due to the added electronics noise that is captured in the Noise Figure metric).  So now the SNR at the output of the amplifier is 17 dB, but the thermal noise is amplified and now much higher (23 dB higher!) than the actual thermal noise floor present at the amplifiers output (we have smothered it with an amplified version of the noise). The attenuator can now push both amplified waveforms (signal and noise components down by it's attenuation or 10 dB. Thus the signal that we increased by 20 dB is now reduced by 10 dB, and the noise that we increased by 23 dB is reduced by 10 dB. However the underlying thermal noise is now closer to the amplified noise level, 13 dB away, contributing a small amount to total noise (explaining the additional 0.2 dB in the 3.2 dB Noise Figure result).
Notice the power summation of the amplified thermal noise with the thermal noise 13 dB lower is a total noise increase of 0.2 dB (it all adds up, demystifying cascaded noise figure!):
$$10 \log_{10}(1+10^{-13/10}) = 0.2 \text{ dB}$$


IIR initial condition / state meaning I am studying filters, in particular the Infinite Impulse Response ones.
I digged into theory and practice and I gain some decent knowledge on the topic. However, I can't really figure out or find good resurces about the so called "Initial conditions" or "State of the IIR".
So far I kind of understood that it is related to the fact that we are using a feedback on previous outputs (y(n)), but I couldn't find a clear and exhaustive explanation on the topic.
In the Digital Signal Processing, Proakis at page 91 of the book (107 if considering the PDF index) it says:

The term y(n_0 -1) is called the initial condition for the systems....

In Digital Signal Processing with Scipy at page 5, the author shows an example of a low pass IIR filter and shows how the initialization of the filter state changes the result of the output.
In the latter paper, the author uses 4 values as initial state, and I can't figure out why.
I'm having some problems in understanding the meaning of these initial conditions /state, where do they come from and how should I compute them. I know that Scipy and Matlab already provide functions to compute those states, but I would like to have a more deep knowledge and understanding of what's going on behind the "scenes".

The states depend on the form of an IIR filter. Take second order filters for simplicity. The following picture shows the structure of a transposed direct form II filter.

The differential equation is
$y(n) = b_0  x(n) + s_1(n)$
and the filter states can be calculated by recurrence
$s_1(n) = b_1 x(n-1) + s_2(n-1) - a_1 y(n-1)$
$s_2(n) = b_2 x(n-1) - a_2 y(n-1)$
As for (transposed) direct form I, it should have a filter states with 4 elements. The following is the transposed direct form I implementation of a second order filter.


Audio Fingerprinting I am doing the project in audio fingerprinting, in that after I convert the signal into mono, I have to apply the hanning window, the size of the window is 370ms, and I have to overlap the signal with a factor of P=28/32.

How can i do that?

What is meant by 28/32?

How can I overlap this signal with this overlapping factor?

How can i do that?

Generate a 370 ms long Hann window, take a 370 ms chunk of your audio, then multiply them together, sample-by-sample.  You're just taking a chunk of audio and fading it in and out at the ends to get rid of any glitches at the start and stop of the chunk.

  What is meant by 28/32?

I'm not sure.  I'd guess it means that the waveforms overlap by 28/32 = 87.5%.  This just means that the chunks of audio you read should overlap with chunks you've already read.

Presumably your sample rate is 44100 Hz, and the window length is 214 = 16384 samples = 371.5 ms?  So

Generate a 16384-point Hann window and store it in memory
Read audio from sample 0 to sample 16383
Multiply them together to get a windowed chunk
Do FFT or whatever other processing on the windowed chunk
Read audio from sample ~2048 to 2048 + 16384 (87.% overlap with previous chunk)
Multiply with window to get a windowed chunk
Do FFT or whatever
Read audio from sample ~4096 to 4096 + 16384
...

The numbers are probably off-by-one, but you get the idea.


DCT - Coefficients and Basis Function I am trying to understand how the Discrete Cosine Transformation works but I am not sure if I am at the right road.

Assuming that I have an 8x8 pixels image and I am applying the DCT to this sample. What I will get as a result is a 8x8 matrix of the DCT Coefficients.

In order to convert from this to the original image I will need to multiply each DCT Coefficient with its corresponding basis function and sum the result. But how do I calculate the basis function if DCT only returns the coefficients?

When DCT is defined by a matrix, then this matrix contains the necessary information to build the basis functions.

Suppose that $I$ is your $8\times8$ block, and $D$ a real $8\times8$  matrix for a 1D DCT (with column-wise vectors). Then $D^TI$ applies the DCT on columns, and $ID$ does it on rows. Thus, a 2D DCT yields a  $8\times8$ matrix $C$ of coefficients defined by:

$$C = D^T I D\,.$$

Consequently, we have $$I = D C  D^T  \,.$$

If we call $\Gamma_{m,n}$ the matrix such that every element is zero, except for $\Gamma[m,n]=1$ (the canonical basis), then

$$C = \sum_{m,n} C[m,n]\Gamma_{m,n}$$

hence, by linearity:

$$I = \sum_{m,n} C[m,n] D \Gamma_{m,n} D^T  \,.$$

Each matrix $D \Gamma_{m,n} D^T $ is an element of the basis you are looking for.


Frequency response of discrete PT1 in MATLAB I want to calculate the frequency response of a discrete PT1 system with T=T1 and K=1. (the code is at the bottom of the question)
To validate the magnitude response, I want to use an sinusoidal input signal with linearly increasing frequency from 0 to fs/2.
The magnitude response should be the envelope of the transient response.
However, this is only the case if I change the line

u = sin(2*pi*f.*t);

to

u = sin(pi*f.*t);

Why? I assume that I have any kind of faulty reasoning here, but I can't find it.

Thanks in advance.

Full code:

clear

%% Input (linear frequency sweep)
Ts = .1e-3;
duration = 30;
t = (0:Ts:duration);
fs = 1/Ts;
f = (fs/2)*t/duration;
u = sin(2*pi*f.*t);

%% PT1-system continuous
T1 = 1e-3;
G = tf(1,[T1 1]);

%% PT1 discretized
H = c2d(G,Ts,'zoh');

%% PT1 trasient response
x0 = 0;
x = lsim(H,u,[],x0,'zoh');

%% PT1 freq. response
[h,~] = freqz(H.Numerator{1:end},H.Denominator{1:end},f,fs);

%% Plot
figure(1),clf(1),hold on
plot(u),plot(x),plot(abs(h))

Found the error by myself now.

I did not use the instantaneous frequency for the calculation of the sweep signal, as stated here or here.

The corrected code is:

clear

%% Input (linear frequency sweep)
Ts = .1e-3;
duration = 30;
t = (0:Ts:duration);
fs = 1/Ts;
f = (fs/2)*t/duration;
finst = 1/2*(fs/2)/duration*t;
u = sin(2*pi*finst.*t);

%% PT1-system continuous
T1 = 1e-3;
G = tf(1,[T1 1]);

%% PT1 discretized
H = c2d(G,Ts,'zoh');
x0 = 0;

%% PT1 trasient response
x = lsim(H,u,[],x0,'zoh');

%% PT1 freq. response
[h,~] = freqz(H.Numerator{1:end},H.Denominator{1:end},f,fs);

%% Plot
figure(1),clf(1),hold on
plot(u),plot(x),plot(abs(h))


What distance metric can I use for comparing images? I usually use the mean squared error (MSE) or peak signal-to-noise ratio (PSNR) to compare two images, but this isn't good enough. I need to find a formula that returns a very big distance between an image A and its pixellated (or blurred) version B, but I don't know how to proceed. What would be a good metric for my needs?

The following is not intended to be an answer, but is a statistic that will help us choose an appropriate image comparison technique based on the characteristics of the images you are analyzing.

The first step is to plot a "delta histogram" as follows:

for (x,y) in [0, width] x [0, height] begin
    delta = abs( SecondImage(x, y) - FirstImage(x, y) )
    hist[delta] += 1
end

Given a plot of this histogram, we will know a bit more about the "magnitude" of the changes you are looking for, and will make better recommendations.

(Alternatively, post some sample images. Remember that if the sample images aren't representative of the image differences you are interested in, we might make inferior recommendations.)

You can also test Structural Similarity (SSIM) on your image set and post your results here. Remember that SSIM is designed to mimic human's ability to recognize obstructiveness of image degradation, so it would detect pixelation but maybe not blurring.

If your images are not photographic images (or, are scientific images that are not ordinary subjects of photography), then please also post examples of their 2D autocorrelation, suitably cropped and scaled.

Face recognition is too big a topic to be discussed in a single question. Blurring arises in multiple context in face recognition - it can be a data quality issue, or it can be done intentionally as an intermediate step in data processing.

In face recognition we want to detect the identity of faces, therefore we have to ignore image differences that are not caused by identity differences. The basic category of differences that should be ignored in face recognition are: pose, illumination and facial expression.

A general approach to ignore irrelevant differences is called normalization, which attempts to apply various operations and transforms on the input image to obtain a "canonical" or "preprocessed" image, which in turn can be used for identification.

A second approach is to extract features from images that are highly-invariant from the irrelevant factors.

The quality of a face image is subject to the capturing device and the environment where it was captured. When a face image is captured without cooperation of the subject (such as from a security camera), poor image quality is an unavoidable consequence and had to be remedied by software so as not to hamper identification.

In cooperative capture, a computerized measure of image quality is good: the operator can be notified of quality problems and the image can be re-taken.

Blurring can also be an example of malicious tampering of biometrics in order to evade detection (along with occlusion and disguise). If the image is encoded digitally, a digital checksum and cryptographic signature is sufficient to solve the problem completely. If the blurred image is submitted in physical print by an impersonator, a computerized measure of facial image quality can be used to reject such submissions.

The lack of 2D-localizable features or interest points in a certain part of facial image can be a sign of intentional blurring.

However, the broad category of digital image tampering (by a skilled user of image-editing software), can only be dealt with Digital image forensics which compares pixel statistics against known camera models.


How to determine the filter type based on the impulse response? When given a discrete impulse response, let us say it is from an FIR filter so it is finite, how can I assume the filter type?
The task is to recognize or to preclude one or more of the following types:

Low pass
High pass
Band stop
Band pass

For example, I read that a zero at $h(n = 0)$ (where h is the impulse response of my filter) means that the filter is a low pass but how can I say that with precluding the band stop immediately?

If you know for sure that the filter type is one of the four given ones, then the exercise is quite straightforward. Look at the zeros of the corresponding transfer function

$$H(z)=\sum_{n=0}^{N-1}h[n]z^{-n}\tag{1}$$

where $N$ is the filter length. The zeros at DC ($z=1$) and at Nyquist ($z=-1$) are given by

$$H(1)=\sum_{n=0}^{N-1}h[n]\tag{2}$$

and

$$H(-1)=\sum_{n=0}^{N-1}h[n](-1)^n\tag{3}$$

It should be obvious that a (classic) low pass has a zero at Nyquist but no zero at DC. For a high pass filter the opposite is the case. A band pass filter has zeros at DC and at Nyquist, whereas a band stop filter has no zeros at DC or at Nyquist.

PS: What you read about the zero of the impulse response is either non-sense, or you misunderstood what they meant.


removal of blackman harris window results in signal higher than expected We are capturing data using a PX14400 adc.  The samples are 16bit unsigned integer, 200M samples/sec, I data only.  100MHz bandwith.  We have a cw signal at 51MHz, -11.4dBm.  The -11.4dBm is measured by a power meter so that includes all losses right before being digitized.

Question 1:  when I apply a blackman harris 4 term window, the fft data calculates to -14dBm which sounds reasonable since bh loss is about 3dB, however when I dont apply the window the fft data calculates to -5.3dBm.  I would have expected a 3dB jump, not almost 9.  Why the large jump?

Question 2:  Does anything jump out in the code below that might be off in my calculations.  For brevity not all the code is shown (threads, board init, etc).

// Initialzied elsewhere in the class
const int fftSize = 8192;         // size of the fft
complex<double> array[fftSize];   // array of complex (r+i)...used post fft
double dData[fftSize];            // array of doubles, used pre and post fft
fftPlan = fftw_plan_dft_r2c_1d(fftSize, dData, reinterpret_cast<fftw_complex*>(array), FFTW_ESTIMATE);

void processNewBuff(uint16_t *rawdata) {

        // convert the samples to V.  This formula is provided by the vendor
        // dData is an array containing the converted to V data (double)
        // rawdata is an array containing the initial adc counts
        // since we only fft the first 8192 values, we effectively throw out the rest of the data
        // we skip the rest of the 10ms sample (200MS/s * .01 = 2MB sample, 8192 is processed the rest is discarded
        double dInputRange = .220;
        for (int x=0; x<fftSize; x++) {
            dData[x] = (-dInputRange)/2 + ((rawdata[x]/65532) * dInputRange);
        }

        // apply blackmanHarris
//        double bhtmp;
//        for (int x=0; x<fftSize; x++) {
//            bhtmp = (double)(2 * 3.1415 * x / fftSize);
//            dData[x] *=(double)(0.35875 - (0.48829*cos(bhtmp)) + (0.14128 * cos(2*bhtmp)) - (0.01168*cos(3*bhtmp)));
//        }

        // do the fft.  fftPlan is a "fftw_plan_dft_r2c_1d", discreet fourier transform, real 2 complex
        // see: http://www.fftw.org/doc/Real_002ddata-DFTs.html
        // after the fft, we throw away the upper half of the fft since we only started with 'I' data
        fftw_execute(fftPlan);

        // normalize the fft for power
        for (unsigned int x = 0; x <(fftSize/2); x++)
            array[x] = array[x] / complex<double>(fftSize, 0);

        // convert the fft to power spectrum
        // I have this separated into 3 distinct lines for readability, final code will be all in a single line
        for (unsigned int x = 0; x <(fftSize/2); x++) {
                dData[x] = sqrt(real( array[x] * conj(array[x]) ));
                dData[x] = ((dData[x]) / 50) * 1000;     // convert to W, then mW
                dData[x] = (20 * log10(dData[x]));
        }
}

The reason for that is that you don't normalize the DFT samples properly. Dividing by number of samples in time domain is valid only for Rectangular Window. For simple case of DFT, you should divide your amplitude by sum of window samples:

$$S=\sum_{i=0}^{N-1}w_i $$

This for Rectangular Window will be obviously equal to $N$ ($8192$ in your case), whereas for Blackman-Harris Window it's about $2939$, therefore the Coherent Gain is about $0.36$ instead of $1$.

For reading more about this interesting topic I recommend you these papers (especially the first one):

  G. Heinzel et al - Spectrum and spectral density estimation by the
  DFT, including a comprehensive list of window functions and some new
  flat-top windows - especially chapter 9

  F. J. Harris - On the Use of Windows for Harmonic Analysis with the
  Discrete Fourier
  Transform


What is imaginary part after IFFT and how to use it? let's say we have signal in time domain.
After we make FFT on it we will receive that signal in frequency domain.
And as far as I know the real numbers mean magnitudes of each frequency bin, and imaginary numbers mean phase shifh of each frequency bin.
And then if we make inverse FFT we get again the the signal in time domain.
And as far as I know that signal is expressed by real result of IFFT.
But then what mean imaginary part of IFFT output?
And - what is also important - how to use it?
I am asking because I have problems with phase synchronising after IFFT so I wonder that maybe that imaginary part of IFFT has something to do?
What is very strange for me is that I have exactly the same algorithm and I get totally different results on two various computers.
I created the clean sine wave oscillator. And I use it to test my FFT->IFFT algorithm.
And on one machine it seems like phase is synchronised properly if my FFT->IFFT size is 4 times buffer size. I mean that for each one buffer (let's say size 512) I calculate FFT->IFFT size 2048. Where first 512 samples are my sine wave, and the rest 1536 samples are just zeros. And I get quite good results.
But on other machine with exact the same algorithm, the same buffer size 512 and the same sample rate 44100 I have phase issue that cause my signal sound awful.
I have no idea why it is happen. The only thing which come to my head is that the machine with awful IFFT signal is just slower and the awful sound is not caused by unsynchronised phase but because of gaps between each buffer. But I can't believe it while both machines are similarly powerfull.
So finally the main question is how to ensure phase is always synchronised after IFFT. And how to do that to work on all machines the same?
For any help great thanks in advance.

And as far as I know the real numbers mean magnitudes of each frequency bin, and imaginary numbers mean phase shifh of each frequency bin.

No. The magnitude of the complex spectrum $|X(k)|$ means magnitude of each freqency bin, and the phase angle of the complex spectrum $\angle X(k)$ represents phase shift.

But then what mean imaginary part of IFFT output?

If the input of IFFT is conjugate symmetric, the output of IFFT must be a real signal, and vice versa, FFT of a real signal must be conjugate symmetric.


confused in practical example of shift invariance? I am reading book of w smith ,"The scientist and engineer guide to digital signal processing" 2 edition
I am trying to understand topic of shift invariance from theory given on pg 92

I have attached screen shot and i have indicated the example with in red marks.

I am confused in whether shift in variance ,author is talking about is in regard to which system?telephone line or digital filter?
What is our target system in this example that needs to be shift invariant?Is it telephone line or digital filter?

They talk about the shift-invariance (or variance) of the digital filter. In general, the response of a shift-invariant system is independent of a shift of the input signal. The corresponding output is just shifted by the same amount as the input. So if $y[n]$ is the response to an input $x[n]$, then the response to a shifted input $x[n-n_0]$ is $y[n-n_0]$. If the independent variable is time, then such systems are also called "time-invariant".

Naturally, an adaptive system cannot be shift (time)-invariant, because it changes its characteristics over time, and consequently, it will react differently to the same input applied at different times.

The example in the text refers to the fact that not all problems can be solved by shift-invariant systems only; in certain situations one needs adaptive systems.


What are two different continuous-time signals when sampled at $10\textrm{ Hz}$ would produce the sequence $x(n)=\cos(n\pi/8)$? The answer is:

  $\cos(1250\pi t)$ and $\cos(21250\pi t)$

I’m not able to verify.

My understanding is that when you periodically increase $t$ by $1/10$ ($10$ here is the frequency, so $1/10$ is time period) from $0$ to the signal $\cos(1250\pi t)$, you should get $x(n)=\cos(n\pi /8)$.

However, $\cos(1250\pi (1/10))$ is not same as $\cos((1)\pi /8)$.

So I would like the verification and what's wrong with my reasoning? Any help is appreciated.

Considering the aliasing, here's the simple derivation: Given a cosine-signal $x(t)$ by

$$x(t)=\cos(2\pi f t)$$

which we sample with sampling frequency $F$, i.e. $t=n/F$, we get

$$x[n]=\cos(2\pi f/F n)$$

and if we now factor $f=kF+f_0$ with $k\in\mathbb{N}$ you'll see that

$$x[n]=\cos(2\pi (kF+f_0)/F n)=\cos(2\pi f_0/F n)$$

For an illustration of aliasing, you can have a look at one of my articles on aliasing.

So, cosines that are an integer multiple apart in frequency, yield the same sampled signal.

Now, to have

$$cos(2\pi f_0/F n)=\cos(\pi n/8)$$

You need $2\pi f_0/F = \pi/8$ and hence $f_0 = F/16$. Hence, your proposed solutions are not correct in my opinion, because $1250\pi t=2\pi 625=2\pi(62*F+5)$.


What's the number of FFTs in Welch-based PSD? I've read from the book "Proakis, J. G., & Manolakis, D. G. (2004). Digital signal processing. PHI Publication: New Delhi, India." Chapter 12 that the number of FFTs required for estimating the PSD from Welch method with 50% overlap is $2N/M$ where $N$ is the amount of data that you have from the signal, $M$ the size of the segments you choose for the method and I guess the $2$ is for the 50% overlap but from the articles I've read they usually do one FFT less than the one described in the equation, why is that? or what am I missing from this interpretation? My theory is that maybe the FFT missing it's supposed to be done with the last 50% of the data and another 50% of data completed with zero padding. Am I correct?

It should really be 1 less than 2N/M. There's nothing to overlap the last "half-transform" (N/M) worth of samples.
You're also right, one could zero-pad that, but then that wouldn't be an Welch estimate.
So, I'll do something I rarely do: Blame Proakis for this :)


Reconstruction of bandpass filtered signal from decimated version of itself I know how to up-sample a discrete, real signal (by an integer factor n) that is band-limited to frequencies between $f_1 = 0 \mathrm{Hz}$ and and $f_2$: Just insert $n-1$ zeros between every original sample (thereby multiplying the sampling rate by $n$) and send the result through a brick wall (sinc) filter to get the interpolated values.

In my DSP book I've read about a generalization of the sampling theorem:
If I have a discrete, real signal that is band-limited to frequencies between $f_1 = f_0  k$ and and $f_2 = f_0  (k + 1)$, with $k$ even, I can just subsample it with a frequency of $2 f_0$ (i.e. decimate) it, and I will retain the complete information (basically the frequencies get shifted down to the baseband between $0 \mathrm{Hz}$ and $f_0$).

My question to you: does the following method let me correctly reconstruct the original signal (in the frequency range between $f_1$ and $f_2$ as given in the previous paragraph) from the decimated signal?

Just insert $k-1$ zeros between every original sample (thereby multiplying the sampling rate by $k$) and send the result through a brick wall filter (bandpass filter with passband $f_1$ … $f_2$) to get the interpolated values.

If this method is not correct, what needs to be changed to get the correct results?

I assume the problem could also be solved by converting to an analytic signal, upsampling in the baseband and then shifting the frequencies up by $f_1$. I want to avoid this because of the additional cpu load this would cause. Please let me know whether you think the method above would work as well, or some other ways to save some processing costs.

Looking forward to your tips,

Flo

Zero stuffing (which is what you propose) will cause a lot of harmonics in the upsampled signal.  They will be proportional to the strength of the original signal.  Your brick wall filter (a low pass will do, no need for a bandpass) will greatly reduce them, but won't get rid of them entirely.

Depending on your application, the result of zero stuffing and filtering may be fine.  If not, you could minimize the harmonics with a couple of methods.

You could duplicate samples (original sample+ (n-1) copies.)

You could interpolate (calculate values along the line between sample und sample+1.)

Duplicating is easy and will cost about as much processor time as zero stuffing.
Interpolating would be smoother but would cost more time as you must calculate individual values.


value of complex multiplier in DFT/FFT I am studying proakis book digital signal processing using matlab 3rd ed

but i am bit confused about calculation of value complex multiplier W

fig is attached. i am not able to understand how the values of W are found/calculated in this fig in red enclosure

how/why we know that W base 4 and raised to the power 1 is equal to -j

Using Euler's formula you can see that

$$W_4=e^{-j2\pi/4}=e^{-j\pi/2}=\cos(\pi/2)-j\sin(\pi/2)=-j$$

With $j^2=-1$ it should be easy to verify that

$$(-j)^0=(-j)^4=1\quad\text{and}\quad(-j)^2=(-j)^6=-1$$

It's important to develop a geometric intuition which allows you to immediately see these equivalences. Multiplication with $j$ corresponds to a rotation by $\pi/2$ in the complex plane, and multiplication with $-j$ is a rotation by $-\pi/2$. So in general for $k\in\mathbb{Z}$ you have

$$(-j)^{4k}=1\\(-j)^{4k+1}=-j\\(-j)^{4k+2}=-1\\(-j)^{4k+3}=j$$


PLL for Phase Demodulation and Carrier Tracking I've been playing around with using a PLL for phase demodulation and it's working pretty well. The PLL locks to the carrier and the error term is the demodulated message signal. However, depending on what the underlying message signal is, it is sometimes necessary to configure the PLL to have a very small loop bandwidth - the effect is to prevent the PLL's internal state to respond too quickly to deviations caused by the message signal.

On the otherhand, it's very useful for the PLL to automatically handle locking to the carrier at an offset frequency and/or track doppler - in this situation, a larger loop bandwidth is essential.

With this in mind - is it a valid approach to have a PLL with a default acquisition mode with a large loop bandwidth, and then once lock is detected, transition to a smaller loop bandwidth to give good phase demodulation performance? One idea I had for this is computing the frequency estimates via the loop's phase (e.g. $f_{loop}[n] = \phi_{loop}[n] - \phi_{loop}[n-1]$) vs the detected input sample phase ($f_{i}[n] = \phi_{in}[n] - \phi_{in}[n-1] $), if these are about the same (after averaging/lpf) the loop would be considered locked and the bandwidth reduced.

Are there any well known/used approaches to take in this problem - (solutions beyond the single PLL are welcomed too).

thanks.

Yes this is very common to have a dynamic loop bandwidth such that during acquisition the loop bandwidth is wider, and then once acquired to tighten it up for better noise performance.

A typical loop will have an error signal determined which is presented to the input of the loop filter. The filtered version of this error signal can be thresholded and used as a lock determination metric. The loop must also be designed such that the loop bandwidth can be transitioned from wide to narrow without disrupting the lock state.

Note that even after acquisition you would like the loop bandwidth to be as wide as possible depending on the dynamics of the link (you mention Doppler, of interest after acquisition would be the maximum possible rate of change of the Doppler offset), but you don't want it so wide that it starts to track out the phase variation inherent in the modulation itself. The phase noise of your local oscillator (LO) plays a big part too in optimizing this loop bandwidth for carrier tracking; if you make the loop bandwidth too tight then your LO and similar jitter sources (ADC clock) start contribution to your noise in detrimental ways.

To this point consider the plot below showing the effect of local oscillator phase noise on a QAM signal and how it interacts with a carrier recovery loop:

Consider two "slices" from the phase noise spectral density for an example oscillator:

The carrier tracking loop is a high pass filter to these noise contributions; the the lower frequency phase noise contribution (which has a significantly higher noise level if not filtered) is significantly effected by the carrier tracking loop bandwidth.

Here is the result of a simulation I had done where you can see how the phase noise contribution is reduced from the rejection provided by the (higher) carrier tracking loop bandwidth.

I wanted to hammer this point home since your question implied that a lower loop bandwidth is a lower noise result and it's actually a minimization problem between LO phase noise/ clock jitter, and carrier loop tracking noise (and signal removal from tracking out the signal).

I go into more detail of this at this post:

Loop bandwidth for symbol timing recovery


Poles of the transfer function in Z Transform Given that:

$H(z)$ has 4 poles maximum.
$H(z)$ has a pole at $z_1=a+bi$

Given that the impulse response $h[n]$ is:

Symmetric: $h[n] = h[-n]$
Real: $\forall$$n$ , $h[n]$$\in$$\mathbb{R}$

How we can conclude that the other poles are the inverse and the complex reflection of $a+bi$ ?

The transfer function $H(z)$ is the $\mathcal{Z}$-transform of the impulse response $h[n]$. If $h[n]$ is real-valued, i.e., if $h[n]=h^*[n]$ we have

$$H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n}=\sum_{n=-\infty}^{\infty}h^*[n]z^{-n}=\left[H(z^*)\right]^*\tag{1}$$

So if $z_1$ is a pole of $H(z)$, then $z_1^*$ must also be a pole.

If furthermore $h[n]$ is symmetric, i.e., if $h[n]=h[-n]$ we have

$$H(z)=\sum_{n=-\infty}^{\infty}h[n]z^{-n}=\sum_{n=-\infty}^{\infty}h[-n]z^{-n}=\sum_{n=-\infty}^{\infty}h[n]z^{n}=H\left(\frac{1}{z}\right)\tag{2}$$

which means that if $z_1$ is a pole, then $1/z_1$ must also be a pole.

In sum, if $h[n]$ is real-valued and symmetric you get for each complex-valued pole $z_1$ three additional poles at $z_1^*$, $1/z_1$, and $1/z_1^*$


Shifting in DFT (FFT) Input image

Listing-1

i = imread('Untitled.png');
i = rgb2gray(i);
F = fft2(i);
%%%F = fftshift(F);
F = abs(F);
F = log(F+1);
F = mat2gray(F);
imshow(F);

Output

.

Listing-2

i = imread('Untitled.png');
i = rgb2gray(i);
F = fft2(i);
F = fftshift(F);
F = abs(F);
F = log(F+1);
F = mat2gray(F);
imshow(F);

Output

Seeing the the above two outputs, can you answer the following questions,

Why does the FFT of an image produce such an spectrum where zero frequencies are at the corner of the image?
Why is that a problem (or, is that)?
Why does the shifting operation fixes that phenomenon?

It is all about defining the support of the DFT.
The DFT is periodic with period of $ 2 \pi $.
Hence on any support of the length $ 2 \pi $ you see all the data.

On default, when using MATLAB's fft() the result is on the range $ \left[ 0, 2 \pi \right] $, namely the Low Frequencies are at the beginning of the output vector.
Using fftshift() the DFT is shifted to match $ \left[ -\pi, \pi \right] $ support which means the DC Frequency is centered.

The same logic works on 2D DFT.


How can i derive step response in terms of impulse response from the convolution sum? If it is discrete LTI system,
How can i derive step response in terms of impulse response from the convolution sum?

In an LTI system, any linear operation on inputs, is directly imposed on the outputs, that is, if an LTI system responses to $\delta[n]$ as $h[n]$, then it responses to $u[n]=\sum_{k=-\infty}^{n}\delta[k]$ as $s[n]=\sum_{k=-\infty}^{n}h[k]$.

This is also clear from the convolution since $$s[n]=u[n]*h[n]=\sum_{k=-\infty}^\infty h[k]u[n-k]=\sum_{k=-\infty}^n h[k]$$


Frequency Shifting in MATLAB I am currently attempting to demodulate an audio file using MATLAB. After taking the fft of the signal, I receive this:

After some more coding, I found the following values:

Fs = 44100;
Window size = 441000 (by taking length(signal))
Fc = 11416

My guess is that this signal is a SSB with a full carrier, and I have read that to recover the original signal from the IF SSB signal, the SSB must be frequency shifted down to its original range of baseband frequencies.

The question here is: How could I shift the frequency spectrum of this signal to center at 0 Hz?

I have been reading up on the fft size/bins and its relationship to the sampling frequency. Will this value be needed? I have also tried using fftshift(), but to no avail. The goal here is to shift the frequencies to be centered at 0 Hz, use a filter to get rid of some of the frequencies, then use an ssb demod function to hopefully get the correct sound from it.

Thanks in advance, I appreciate you taking the time to help another guy out.

Try to remember what happens when you multiply a signal with a cosine function with a frequency of Fc.

Namely, what's the shifting operator in Frequency domain?

If we're talking on visualization, so the operator in MATLAB is called fftshift()

Enjoy...


Confirming what a compression ratio is I have looked at various sources, but still have not been able to answer the question: what is compression ratio? I have come to the understanding that this is the definitive answer:

compression ratio = (block size * bits per pixel) / (dct coefficients * bits per pixel)

--or--

compression ratio = amount of bits in original / amount of dct coefficients

Please confirm.

Neither.

A true compression ratio is: "original file size in bits" divided by  "compressed file size in bits". A practical (based on disk limits) compression may embed the chunk size effect: "original file size in number of chunks" divided by  "compressed file size in chunks", less favorable.

Some of the main reasons for "neither" are:

DCT is not the only method for image compression: dct co-efficientswould be meaningless for non-DCT methods
in JPEG, DC coefficients are further predicted in a raster fashion
images have color components, further decorrelated and differentially quantized to reduce variability
even assuming that DCT can stand for "transformed into less redundant symbols": you are loosing the header, the entropy or dictionary coding, and never complete the construction of an actual string of compressed bits.
the "amount" of DCT coefficients, or "dct co-efficients" is not precise enough to be counted.

However, an estimate of the number of bits taken  by quantized DCT coefficients (through entropy for instance), when divided from the original image size, can provide a rough estimate of what you can expect from JPEG compressing, at least at intermediate bit-rates.


Sampling frequency of audio Based on what we choose sampling frequency when sampling some song? Do we use Nyquist–Shannon sampling theorem, or maximum sampling frequencies are pre-defined, so it will never come to aliasing effect?

Of course you have to obey Shannon-Nyquist when sampling anything.

But, of course, sound cards were designed by engineers, so they have analogue filters that avoid aliasing at the sampling rates the cards offer.

So, if you're using a sound card, any sampling rate it supports should be largely alias-free. Use the sampling rate sufficient for your signal of interest and comfortable for your application and available from the sound card hardware. In effect, that means you can chose between 44.1, 48 and 96 kHz on most sound cards.


what is the periodicity of $e^{j\omega n}$ when $\omega=0$? what is the periodicity of $e^{j\omega n}$ when $\omega=0$? I think it is $N=1$. But signals and systems book by Oppenheim and Willsky page no.28, mentions it as undefined.  I can understand this for continuous time signals but for discrete time case, I think unconditionally it is $N=1$ because $x(n) = 1$ and so $x(n) = x(n+1)$ for all $n$. What do you think?

Mathematically speaking, the constant function is periodic with any period $P>0$, because is satisfies $f(x)=f(x+P)$ for all $x$.

The existence of a fundamental period depends on whether $x$ is discrete or continuous. In the first case the smallest positive period is $P=1$. In the second case there is no smallest period and therefore the function does not have a fundamental period.

In any case, claiming it is "undefined" is not accurate. This is all very well defined.


What is meant by Two pass FIR Filter?? (A basic question) I have a confusion about what does a two pass FIR (bandpass) filter with order 40 means?? Passband frequencies are [8 13]. Type 2 FIR is same as two pass??

I have check some previous literature which shows Linear-phase FIR filter can be divided into four basic types.

TYPE I     symmetric length is odd

TYPE II    symmetric length is even

TYPE III   anti-symmetric length is odd

TYPE IV    anti-symmetric length is even

But according to mathworks ([firtype][1])   [1]: https://jp.mathworks.com/help/signal/ref/firtype.html

Type 1 - Symmetric factor with even degree

Type 2 - Symmetric factor with odd order

Type 3 - Asymmetric coefficient with even order

Type 4 - Asymmetric coefficient with odd order

kindly correct me if my understanding is wrong : Length = (Order-1)/ 2

If Type II is two-pass then the order of the filter will be even or Odd??

You can't in general categorize a two-pass filter in any of those listed filter types. "Two-pass" means that the filter processes the data in two passes; the output of the first pass is used as the input of the second pass.

The two-pass filter's impulse response is the convolution of the impulse responses of the filters in the two passes. This comes from the algebraic properties of convolution:

$$\text{output}= f*(f*\text{input}) = (f*f)*\text{input},$$

where $*$ denotes convolution and $f$ is the impulse response of a single pass. The two passes have an identical impulse response to match your description. Because convolution is equivalent to multiplication in the frequency domain, the two-pass filter's frequency response will be the square of the frequency response of a single pass.

With infinite impulse response (IIR) filters, also a configuration where the second filter processes the data backwards can be useful, and the two filters may be connected in parallel so that the outputs from the two passes are summed to form the composite filter output. Such configurations are not useful with FIR filters.


Avoiding identifying objects within previously identified objects in MATLAB Apologies if the question title is confusing.

I've got some code below in MATLAB, which identifies squares/rectangles in an image.

img = imread('image.jpg');
imgbw = im2bw(image);
filter = bwareafilt(imagebw, [800 10000]);
fill = imfill(imgbw,'holes');
area = bwareaopen(fill,15);
final = bwlabel(area);
dets = regionprops(final,'boundingbox');
figure(1), imagesc(I1); hold on;
for count = 1 : numel(dets)
    gg = dets(count).BoundingBox;
    rectangle('position',gg,'edgecolor','r','linewidth',2);
end

I keep getting results where I am identifying squares, within another square I have already identified. i.e.:

What I wanted, was to identify the yellow as one square/rectangle, and leave it at that.  But as you can see, I've identified my rectangle/square, but also identified squares in the middle.  Does anyone know how I can get ignore shapes, that are also inside previously identified squares?

There are two alternatives to this.

Adjust your opening parameter so that you end up with a fully connected shape. Say, instead of 15 you can use 25.
You can find the areas of all intersecting polygons using e.g. this. If the area of an intersection is higher than a given threshold, discard the polygon.

After finding overlapping polygons, you could consider averaging the hypothesis so that you end up with a more correct shape pose.


Good methods for hardware based 1-2 word speaker independent speech recognition? I want to detect, in real-time, a 1-2 word spoken phrase within continuous speech, speaker independently, in hardware. The smaller the number of gates the better.

What is the best way to do this? Are the requirements feasible?

There is really no need to do this in hardware. You can do the same on any cpu and it will be more flexible, energy efficient and updatable. You can pack simple computing core into your chip and run on it if you want. I don't think you'll gain anything from hardware implementation.

For a good implementation of keyword spotting you can check CMUSphinx http://cmusphinx.sourceforge.net, it has recently introduced good efficient keyword spotting. Fixed point arithmetics, noise suppression, tunable threshold. All that you need in single package.


Does 1 kHz sine tone means $\sin(2(1000)\pi t)$ or $\sin(2(500)\pi t)$? Does 1 kHz sine tone means $\sin(2(1000)\pi t)$ or $\sin(2(500)\pi t)$?

The trigonometric functions "do not know" what a Hertz is and they do not care either. The only thing they know is that a full circle is $2 \pi$ radians. Whether this circle concludes in days, hours, picoseconds or a slice of it represents the angle a force is applied to some lever, is immaterial.

$2 \pi \omega$ expressed in Hertz, denotes a rate. A rate of going around a circle at the time span of a second. $y = \cos(2 \pi 1 t)$ where $t$ is in seconds, would have concluded 1 circle, composed of $2 \pi$ radians, by the time $t$ ticks to 1.

To make it conclude the circle faster, we multiply the "passing of time" (denoted by $t$) by some number $f$.

Therefore, a 1kHz tone is $2 \pi 1000$ radians per second.

Hope this helps.


How can I smoothly interpolate between 2 position? I've got a 1D signal (position of a servo motor over time) and I've extracted 'peaks'/'key' positions picking running average "local extrema" points.

Below is are 2 plots from 2 servos and the white markers indicate key positions I'd like to interpolate between:

I've played with linear interpolation, but the motion looks robotic. I'd like to get a more natural/smoother motion and I imagine it's posible by using bezier interpolation between two key points. These two points would be the start and end point of the bezier I imagine, but I need to solve/find out the 'anchor' points as well. Since I have all the points I can get the velocities as well I imagine, but not sure how to plug them into the equation.

Any idea on how I can find the other two members of a bezier equation for easing/smooth interpolation ?

Also, since I'm new to this: is my approach ok ?  Is there a simpler solution I'm missing ?

When you talk about bezier curves, it sounds like you think about them from an "illustrator" point of view, which is not totally right when it comes to spline interpolation (most probably what you are looking for).

Splines are piecewise curves that pass thrhough points. Bezier curves are third degrees splines between two points, and a series of them can form what "looks like" a single curve in a drawing aplication, but are actually lots of individual curves tied together (it counts as if the intermediate knots are doubled).

Illustration concepts aside, I think you could choose any 1-D spline interpolator from some popular library. I use Python, so I'd go with Scipy, take a look:

http://docs.scipy.org/doc/scipy/reference/interpolate.html#d-splines


Simple $\sin(2\pi 1000t)$ Fourier transform in PSpice not behaving as expected So I have this very basic circuit show below which I am simulating with PSpice.

Now, when doing the Fourier transform of $\cos(2\omega1000t)$, I expect to see two impulses(one at -1kHz and one at 1kHz), at half the amplitude of the cosine. The simulated results, however, shows the impulse at 1kHz with an amplitude of just below 1V. This is shown below:

My question is, why is this the case? Why is the amplitude not 500mV? I'm probably just missing something trivial.

Thanks in advance!

(Side note - why is there no pspice tag? Is this the correct form to post this question?)

PSPICE is showing the "Double Sided" Spectrum with the assumption (as it is) of a real signal. This means the spectrum is positive frequency axis only and the result will be double of what you would expect to see with a "Single Sided" Spectrum, basically matching the amplitude of the sinusoidal component represented.


square wave frequency representation New to this signals stuff and i'm confused about the frequency representation of the square wave. Correct me if i'm wrong, a periodic square wave is composed of odd harmonics sine waves which are multiples of the fundamental frequency, In the frequency domain this just represents decreasing impulses at every odd multiple of the fundamental frequency . What relation does the sinc function have with the square wave? and which is the correct spectral representation of the square wave? Thanks!

The transform of a single rectangular pulse is a Sinc function.  Each added translation of that pulse adds a twisted repitition of that Sinc.  The twist will equal 2 pi at all spectral distances that are odd integer multiples of the reciprocal of the repeat period.  Repeat, and eventually all the twisted portions will all cancel out (form circles in the complex plane around zero) except at those points representing the odd harmonics.  Thus an infinite repeating square wave will be left with only non-zero spectrum at odd harmonics in its transform.


Using MATLAB's Spectrogram Function for Analysis In the documentation - sepctrogram().
Question 1)s = spectrogram(signal) and spectrogram(signal) are two commands to plot the spectrogram. However, the variable s is complex valued. I am unable to understand which output of the spectrogram is used to generate the image plot?
Question 2) How to determine the best values of the parameters window and noverlap ? Should noverlap be 50% of the signal length (number of elements in the time series) or 90% etc? What if it is zero then what does it mean? My dataset has sampling time = 1sec. I remember reading somewhere that the window should be at least roughly twice as long as the period of the lowest frequency. So, for my case is w=2 since frequency = 1?
I was thinking of using pspectrum(signal,'spectrogram') which outputs the spectrogram and use the output values as inputs to the spectrogram() function. But again, I don't know which output values from pspectrum can be used, if at all that is possible.

See the MATLAB documentation: s = spectrogram(x) returns the short-time Fourier transform of the input signal, x. Each column of s contains an estimate of the short-term, time-localized frequency content of x.. Namely each column of the matrix s is the result of an fft() on some samples of the input. So the plot you see is the magnitude of the columns of s.

Spectrogram is about analysis of Non Stationary signals. So something is changing over time which means it makes no sense to look on the DFT of all samples. The window length is the time you think the signal has the same properties over time. The overlap time should be similar to the time the signal is changing. Something like the time of Fade Out / Fade In, the transient length.

Example
The following code will recreate the figure form the function (Up to the Colorbar and the units of the Axis):
t = 0:0.001:2;
x = chirp(t,100,1,200,'quadratic');

figure();
spectrogram(x,128,120,128,1e3);

s = spectrogram(x,128,120,128,1e3);

figure();
hA = axes();
imagesc(20 * log10(abs(s).'));
set(hA, 'YDir', 'normal');
```