Scale and Rotation invariant color based feature descriptors SIFT/SURF/ORB etc are grayscale based, if the two grayscale bitmap is very similar, but the RGB channel is quite different, how to classify them? So, I think there should be some color based feature descriptors, also have scale and rotation invariant, like SIFT/SURF/ORB, but for RGB or HSV.

But util now, I cannot found such a feature descriptor ......

There is absolutely no difference in the way these algorithms operate over grayscale and RGB images and colour SIFT algorithms simply take into account the additional information contained in colour.

These algorithms operate on typical images that are represented as matrices of (usually) integer numbers. Each number represents the intensity, or brightness, of light coming in to the imaging equipment from a particular direction. Colour information is (usually) represented with more than one information channels due to the way light works. These additional channels, also look like matrices. The raw form of a typical RGB optical photography image is three grayscale images packed together, OR, one image with 24 bit colour depth (8 bit for every colour channel).

There are however different ways by which you can visualise the RGB colour space. One of them is the Hue Saturation Value that returns one single number representing colour (Hue) another one representing the strength of colour (Saturation) and another one representing the brightness (Value).

The point in colour SIFT is to diminish the effect that colour has on tracking features but do it in the same way across all different colours (colour invariant). This can be done, either by pre-processing the input image (e.g. transforming it to a different colour space) so that the effect of colour is diminished and THEN sending it to a typical SIFT, or by incorporating colour invariant features to the standard keypoint selection process employed by SIFT (thus, cSIFT).

Besides SIFT, colour invariant metrics are also required in image indexing where they have been in use for longer, so you might get some additional information from "colour invariant metrics in content based image retrieval" (or this) type of searches.

Hope this helps.


FFT with fewer outputs than inputs Is it possible to write a FFT function, where there is fewer outputs than inputs? What in mean is, we would start computing 16-point FFT, and our result would end up be 12-point.
If it is possible than do we loose any data in the procoess of reduction?

The FFT is just an implementation of the DFT, which is an invertible operation on any vector composed of complex numbers.
Since invertible operations need to be surjective:
No, an FFT can't have less outputs than inputs.
You are, probably, looking for something that "condenses" your data, e.g. looking for features, or projecting onto a smaller vector space that is sufficient to describe your signals of interest; keywords you might be looking for would be "feature extraction" or "compressive sensing".
However, I think it might be sensible for you to post a new question, describing exactly what you actually want to do (and not focused on employing the FFT itself).


Count the vehicles flow that pass over one bridge I live in an apartment near a bridge. From my window I can see the cars traveling over the bridge as you can see below:

Is there a relatively simple manner to count each car pass over the bridge? Using an android cell phone pointed to the bridge behind a lunette and using a API for image processing like Open CV?

It is relatively easy to create the "sensing element" for the cars but difficult to get a short term valid count for the number of vehicles.

As far as the sensing element is concerned, this is a simple pulse counter. The signal that creates the "pulse", comes from the optical flow of the vehicles as they pass by a region of interest in front of the camera.

You basically need to differentiate between successive frames, rectify the differentiated signal and integrate it with respect to time with the possible addition of a threshold.

When you differentiate between two successive frames (for instance $I_n, I_{n-1}$, where $I$ is a two-dimensional array of pixel brightness values at some color depth and $n$ is the frame number (integer))  you are basically removing everything that is not moving and highlighting parts of the image that appear to have changed between time instances $n$ and $n-1$.

At this point, you don't really care about the direction of the change and this is why you need to take the absolute difference: $Q_n = |I_n-I_{n-1}|$

If you now focus on a small region of interest of the $Q$ series of images, you will notice that when nothing happens, that region is pitch black and when an object MOVES in front of it the region "lights up" because of all the differences. Therefore, what you can do now is obtain the mean value of all the pixels of this small region of interest (over the $Q$ sequence) and observe it in time.

This will be producing a "pulse" every time that something disturbs the region of interest. The duration of the pulse (in time) will be equal to the duration of the disturbance.

Now, if we were to think in terms of "counting cars over a bridge", you can position your region of interest around a light pole. Past the first frame, the actual scene will not matter anymore because the light frame doesn't move (actually it does move because the whole bridge vibrates slightly and the light poles are tuning in but let's assume that this is negligible in your case). If you do this then you can count disturbances under the light pole during the day and also during the night in a similar way (rather than having to rely on the irregular shapes of the car lights).

Fine, so now we have a way to transform passing vehicles into "pulses". What is the problem? Count the pulses and you get the number of vehicles.

Not exactly, so here is where things become difficult:

There are more than one lanes per direction. Therefore, it is quite possible that two cars fit through the region of interest at the same time thus creating just one pulse. Even worse, a car can be masked by a bus (be behind it) and never be counted. The only way to counteract this is to repeat the same ultra simple "counting" process at multiple points on the bridge (multiple light poles) and obtain an average count. (Of course, a car can travel at the same speed as a bus and remain hidden from the perspective of the camera for as long as it remains on the bridge...but there is no remedy for this....from your point of view).
As you can see, this pulse counting is actually counting disturbances, not necessarily vehicles. Here are a few disturbances that have nothing to do with the vehicles: Atmospheric disturbance (smoke, sea evaporation, anything that causes a slight flickering), movement of the camera, birds and light signals. Imagine for example what happens when the bridge lights will come on. Your pulse counter will register a pulse (because before, the area under the light pole was darker and now, it is lighter) but this pulse is going to be extremely brief. The only way to counteract spurious movement somewhere across the image is by the solution to #1 (multiple counting points + average) and the only way to counteract invalid pulses is by setting thresholds. So, very short pulses (in time) are not counted and very weak pulses (in terms of magnitude) are not counted too.

This is the principle. You can do this with openCV, python, (openCV through Python), MATLAB, even off-line video processing.

Hope this helps.


How is wavelet center frequency computed? PyWavelets (1) takes index of max  DFT magnitude, (2) adds 1 to it, (3) divides by domain, which is the range of input values to the wavelet ("support"). I don't see where steps 2 and 3 come from, or why 1 doesn't suffice.
MATLAB's docs suggest only step 1, also here - but closed source, can't verify. If "center frequency" is frequency that correlates the strongest with the wavelet at some reference scale, that's by definition what the DFT finds, so why the extra steps?

Short version: DFT's bin indices are input length-dependent; "center frequency" is measured relative to the function generating the wavelet. Generated length can vary, so must be accounted for to yield 'correct' $f_c$. This requries injecting information from knowledge of the function, as DFT is blind to "absolute duration"; PyWavelets does this via domain, the range of input values to the function. domain=16, so DFT's bin indices are x16 shifted; undo by dividing.

The answer's in DFT. Take Real Morlet, with the $\sigma >=5$ approximation:
$$
 \psi(t) = e^{-t^2/2}\cos(\sigma t) \tag{1}
$$
Suppose we evaluate it for $\sigma=5$, $t\in[-8,8]$; what's it look like?

So what is its "center frequency", or frequency of a sinusoid with which the wavelet correlates most strongly? Try tweaking $\beta$ here. Spoiler: sigma; the exponential is almost ignored, just slightly shifting optima locations (but strongest correlate still is $\sigma$):

As there's no $2 \pi$ in $(1)$, 5 is the angular frequency, so $f=5 / 2\pi=0.796$, or $f^{-1}=1.26$; these match the red cosine, which cycles once every 1.26 units, or 0.796 times every 1 unit.
But in the discrete realm, we've no notion of "time" unless explicitly labeled; all we have is samples. What happens if we take the DFT of red? Do we peak at $5$, or $5 / 2 \pi$? Neither; both of these numbers are "per second", and DFT doesn't understand that. All it knows is that a cosine completing one full cycle from left to right ... completes one full cycle from left to right. It calls this $k=1$.
Then we seek $k$ such that the cosine matches most closely with red. Since red makes 0.796 cycles every 1 unit, and there are 16 units, then it'll make 0.796 * 16 = (1/f) * domain = k cycles in total. That's the index of max(abs(DFT)); we just swap f and k in the equation, inputting k to estimate f.
All said of red above also applies to blue; their peaks match. Putting it all together:
f_center = (t[-1] - t[0]) / k_max = domain / k_max

But this looks different from PyWavelets
Because PyWavelets looked ahead. First, we never expect our wavelet to be flat, so we eliminate the possibility of k_max=0 in case the wavelet averages nonzero:
argmax(abs(fft(psi)[1:]))  # _k_max - 1

But this shifts all k, such that original k=1 is now k=0; so we do +1:
index = argmax(abs(fft(psi_pad))[1:])) + 1  # _k_max

What if $k > N/2$, where N == len(psi)? Negative frequency. This is fine, except we seek the counts of cycles, which are non-negative, so we just find its positive. Suppose N=101; then, there is no Nyquist frequency, only max positive and negative frequencies (equal) at k=50, 51.

_k_max = 100. This is negative k=1. To obtain positive: N - _k_max = 101 - 100 = 1.
_k_max = 50. This is the positive max; do nothing.
_k_max = 51. This is the negative max; make it (+): N - _k_max = 101 - 51 = 50.

This explains
if index > len(psi) / 2:
    index = len(psi) - index  # now we get actual `k_max`

And that's all. Except it's not; PyWavelets does domain / (index - 1), and +2 instead of +1, and -index + 2. The behavior is identical for most, but not all, possible len(psi) and _k_max; grid search. Consider (len_psi, _k_max) = (101, 50). pywt finds k_max=51, which is a negative frequency.

How accurate is it?
Not even within two significant figures. Authors may have mistaken higher precision for higher precision, when in reality k_max gets scaled down by the exact amount psi gets scaled up. (Exception for precision < 5, for which estimation greatly degrades):

Can we do better? -- Yes; above right is the result. See here. Opened a PR, also on k_max.

More rigorously, but briefly: center frequency isn't necessarily at max frequency component.
From Mallat 4.1, if $\psi$ is a wavelet, then $|\psi|^2$ is interpreted as a probability distribution per
$$
\int_{-\infty}^{\infty} |\psi (t)|^2 dt = 1
$$
centered in time at
$$
u = \int_{-\infty}^{\infty} t |\psi (t)|^2 dt
$$
(i.e. expected value, $\mathbb{E}[|\psi(t)|^2]$) and in frequency ("center frequency") at
$$
\omega_c = \frac{1}{2\pi} \int_{-\infty}^{\infty} \omega |\hat{\psi} (\omega)|^2 d\omega
$$
Corresponding variance relations are given, from which we define Heisenberg boxes for time-frequency resolution:

If a wavelet's FT happens to be symmetric about the peak (as with most used wavelets), then center is at peak's location.

Addendum: 'relative to what' is $f_c$? How to interpret?
Relative to "the function", surely, but how to interpret this in terms of signal data we work with? Like with DFT, to find "physical" values, we need to know the sampling frequency, $f_s$. Physical $f_c$ has units of Hertz, so if we know the number of samples that comprise one second, we can relate "frequency relative to samples" (DFT) to "frequency relative to seconds" (more here). Thus:
$$
f_c \left[ \frac{\text{cycles}}{\text{second}} \right] = \left( f_c \left[ \frac{\text{cycles}}{\text{samples}} \right] \right) \cdot \left( f_s \left[ \frac{\text{samples}}{\text{second}} \right] \right)
$$


How to describe correlated noise after the signal is oversampled? The Gaussian noise in discrete signal models is usually assumed to be independent and identically distributed variables (i.i.d.). Does this mean that the signal must be sampled with the sampling rate of 2W (where W is the bandwidth)? If the sampling rate is faster than 2W, the noise becomes correlated. Now I want to generate correlated noise samples in Matlab, but I am quite confused about how to describe the correlated noise theoretically? I also don't know how to generate correlated noise samples in Matlab. Could anyone explain this to me or give me any tips?

Assume that the input to a simple up-sampler (say by $L$) is $x[n]$, hence
$$y[n] = \begin{cases} { x[n/L] ,\qquad n=kL \\ 0 , \qquad \qquad \ \text{  else } }\end{cases}$$
If you compute the auto-correlation sequence, defined as
$$\phi_{xx}(i,j) = \mathbb{E}(x[i]x[i+j])$$
for an iid white Gaussian process, we have that
\begin{align}
\phi_{xx}(i,j) &= \begin{cases} \sigma^2, \qquad & \text{if } j=0 \\
0, \qquad & \text{else } \end{cases} \\
\phi_{yy}(i,j) &= \begin{cases} \sigma^2, \qquad & \text{if } j=0 \text{ and }i = kL \\
0, \qquad & \text{else } \end{cases}
\end{align}
Unlike i.i.d Gaussian, where $\phi_{xx}(i,j)$ is determined by the lag $j-i$; $\phi_{yy}(i,j)$ is not fully determined by the lag. That said, $y[n]$ is not WSS, whereas $x[n] $ is. With respect to correlation, the quantities $\mathbb{E}(x[n]x[m]) = 0$ for $n \neq m$ is due to the uncorrelation assumption. You can say the same thing on $y[n]$, namely $\mathbb{E}(y[n]y[m]) = 0$ for $n \neq m$.


Possible RF Mixing for Increasing Frequency of a Pulse Train I am planning to build a low-power DC/AC inverter circuit along with MOSFETs and the full-bridge configuration. The MOSFETs are going to be controlled via two 4047s on the both side of the inverter.
According to my needs, the switching frequency should be 382 kHz. However, according to the 4047’s datasheet from Texas Instruments; in order the 4047 to work properly, the maximum output frequency is calculated as nearly 270 kHz and it is way below the targeted switching frequency.
Consequently, I’ve thought about mixing the 270 kHz rectangular pulse train with a local oscillator. According to my Fourier transform calculations, frequency of the output signal of the mixer becomes 382 kHz by using 112 kHz oscillator.
Is this solution is applicable?

It will not work. You can't increase the frequency of a pulse train by mixing it with a sinusoid. This figure shows what happens. The spectrum on top is the pulse train at 270 kHz; the spectrum below is the result of multiplying it by a sinuosid of 112 kHz. The purple lines correspond to shifting the original spectrum to the left, and the orange lines are the original spectrum shifted to the right.
As you can see, the resulting spectrum is not a rectangular pulse train.

If what you need is a 382 kHz sinusoid, you will need to bandpass filter the mixed signal. The BPF will have to be pretty narrow with a short transition band.


determine background noise level via cross correlation The setup:
Sound is being played by the speaker on the iphone. At the same moment i record sound via the microphone. Obviously the speaker sound will be recorded in addition to all the external noise.

My goal:
To determine the level of all external sound which is recorded in addition to the sound from the speakers.

(no implementation yet)
So far my approach is to cross correlate the input stream with the output stream to get the delay between both of them and an integer representing the level of similarity at this point, which will be used to estimate external noise. (if the int drops below a given value the noise is too loud)

Has anyone experience with stuff like that, if this approach will be accurate enough to make a guess if the external noise reaches a certain level of disturbance?
Or will i have t do further analysation in the frequency domain via fft for example.

I'm completely new to this kind of problem so be kind with me, if I'm totally wrong ;)

If you want to get an estimate of the noise, do the cross-correlation, and then multiply your transmitted sound by your cross-correlation peak and then subtract it from your received sounds.  Make sure, of course, that you use the cross-correlation to line up the transmitted sound with when it was received.

What remains after you subtract is the noise with some amount of (hopefully small) error.


Frequency analysis to determine low-pass cut off frequency I collected some data for a practical application, where the signal represents force data obtained from an impact of a punch against a force plate attached to a quasi-rigid rig (it moves once the impact occurs). I have a few questions to understand how to deal with my data.

Features of the signal:

the signal (sampling 1000Hz) has a length of 2 seconds, but the impact last about 60ms. The recording is 2s for practical reasons though. Therefore, before the impact a baseline is recorded (values close to zero), and after the impact vibrations occur until data set back to baseline levels. I need to figure out a reasonable cut-off frequency to use combined with a 4th order butterworth filter to properly filter the signal.

When exploring (FFT) the frequency of that signal should I use the whole 5s long signal or only the 30ms of interest? Due to sampling frequency, data points over 30ms maybe not enough?
Vibrations occurring due to the impact are not of interest, but the frequency is likely to be in part similar to the frequency of the portion of the signal of interest. How to deal with this? I could select from before the impact (arbitrary) to when the signal becomes zero after the impact. Then use from that point onwards to explore frequency of the unwanted signal only?
Impacts are at least 10 per subject, when exploring the frequency of that signal, should I average all trials and all subjects out to obtain an average frequency content or should I perform the analysis for each trial and subject indepedently?

Attached a figure of a typical signal (x = time; y = Newtons)

I think you will get almost the same amount of frequency content for both the cases (fft of 60 ms and fft of 2 sec). The time series indicates that only from 0.6 to 1.2 sec there is vibration. But it's better to use the full signal while doing the FFT as there is no fear of losing any data. If you cut some portion from the time series, you have to be careful about choosing the window (use cosine tapering not rectangular window) as edging effect may be introduced into your frequency content.
Sorry, I didn't get the question. Could you please rephrase it?
I think that for all the trials, the time series would be different. Therefore, it is not recommended to average all the time series. It is better to analyze all the time sequences and obtain the frequency contents. If dominant frequencies for all the trials are similar, then the average frequency can be calculated. If not, you can try to find out the reason e.g. inconsistency in the sample material etc.


FFT plot in Matlab, accuracy I have a small input signal of 60Hz sine wave from signal generator, which is corrupted with 50Hz mains supply frequency. I want to measure the amplitude of the 60Hz signal using FFT because it is very small to see in the oscilloscope.
The Matlab FFT code:

y = data;
Fs = 2048;
[r, L] = size(y);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1)))

But the FFT plot doesn't give a sharp peak at 50 and 60Hz. The plot looks like this: . The consecutive points have high and low amplitude alternatively which gives a saw-tooth like plot. Why is it so? Is the amplitude of 60Hz affected by this?

The phenomenon you see is called "window effect" or also "spectral leakage". It occurs because ideally you would like to analyze an infinitely long signal (which would give you ideal peaks in the frequency domain), but actually you analyze a truncated (windowed) signal. In fact you multiply your signal with a window function - in your case a simple rectangular window - which in the frequency domain corresponds to a convolution with the frequency response of the window. This gives you the effect that you see. You can try different window functions, but you will not get rid of the basic effect.

BTW, for the plot you could try using values in dB, i.e. use $20\log (|Y|)$ instead of a linear plot; this will make things look a bit better.


Should I ever pick the continuous wavelet transform over the discrete one? DWT vs CWT vs STFT Of course I mean in terms of numerical algorithms.

I'm reading various papers on the wavelet transform and why it's better than the short-time Fourier one. The reason that's cited more often is that the WT detects transients of known shape better than the STFT.
In other sources I read that the DWT is essentially a faster version of the CWT. But so far all the MATLAB algorithms I've seen for scaleograms and transient detection have been using the CWT. Is there any specific reason for this?

In general, when should I code a CWT and when - a DWT?

Thanks in advance.

Standard 2-band discrete wavelet transforms have some subsampling, thus they tend to be much less precise in localization. Plus, their shapes are limited, and the ones with finite support are slightly dissymmetric (all but Haar), which is problematic for symmetric pulse/edge detection.

If you don't care about computations, and are interested in visualization, matching, feature detection, CWT might be easier.
If you are interested in  denoising, compression, restoration, DWT are often more appropriate.

Yet, there are many useful schemes bridging the gap between DWT and CWT that can be used in dedicated applications.

A reference paper: Paul S. Addison, 2005, Physiol. Meas. Wavelet transforms and the ECG: a review
References for fast CWT implementations
Matlab code: continuous wavelet transform and inverse
C code: CWT Library: Continuous wavelet transform library
C++ code: Wavelet - A library for online estimation of the Continuous Wavelet Transform


Why are we always interested in mean-squared distortion? When dealing with quantizers, and in many other communications problems, the interest is usually on the mean-squared distortion or mean-squared error, rather than mean absolute error or anything else. Gallager explains the reason here, probably perfectly if I were to extrapolate from the quality of the rest of his lectures, but I don't get it.

He says that the mean-squared distortion of the quantizer maps nicely to the mean-squared distortion of the waveform we sample&quantize, and this mapping wouldn't be as nice if we used another metric such as mean absolute distortion. Why is that? Shouldn't the mean-absolute difference between the original and the sampled-quantized waveform also map to the mean-absolute difference of the quantizer? What makes mean-square different in this sense?

Consider a band-limited function $u(t)$ sampled at an appropriate sampling rate $f_s=1/T$ such that it is perfectly represented by its samples $u(kT)$. If those samples are quantized resulting in values $v(kT)$ then it can be shown that the mean-squared error

$$\epsilon=\sum_k\big|u(kT)-v(kT)\big|^2$$

is equal to the mean-squared error between the waveforms

$$\epsilon=\frac{1}{T}\int_{-\infty}^{\infty}\big|u(t)-v(t)\big|^2dt$$

where $v(t)$ is the waveform reconstructed from the quantized samples $v(kT)$. So quantizing in such a way that the mean-squared error between the original and the quantized samples is minimized, guarantees that also the mean-squared error between the corresponding continuous waveforms is minimized. This is not the case with other error measures.


Fast method to decompress huffman and arithmetic encoded data By nature, when decompressing huffman (and I assume also Arithmetic encoded data) we have to read it in bit by bit and then see how it traverses a huffman code tree that has the actual symbol value at each leaf node. Since processors are designed for byte or multiple byte read/write, this seems to be a very slow and inefficient process.

Since huffman coding and many other methods are used to encode image and video data which produces non-byte boundary aligned results, how do we speed up the decompression process? I would assume that instead of reading 1 bit at a time we read several bits or even bytes and then using some clever method, find all symbols in it at once. I think that this is how all codecs for video and image would work to decompress the compressed image data fast. However, this is only my assumption and I do not know of any details. Could someone kindly reveal the secret to me?

For actual implementations take a look at (C) programming codes for decoders.

Basically the process is simple, you create an abstraction such as

bits = get_bits(n)

Where, n is the n bits you want to read. the var bits is left aligned. It extracts it from some byte_buffer - say a 32 bit or 64 bit symbols.

  NOTE: All, practical codecs put major symbols at 'Byte aligned' even though >internal symbols can be bit by bit.

Now for the implementation efficiency part, basically, irrespective of what huffman code will come out to be, (and hence whatever be the required n).

this is how simple get_bits can be implemented:

// At start
frame_buffer = 0x4444
bits_delivered = 0

n = 1
mask = 0x8000
bits = (mask & frame_buffer) >> 32-n
mask = mask >> n
bits_delivered += n
if(bits_delivered == 32)
    { load_framebuffer( &frame_buffer) }

(Actual code needs to take care of few more stuff boundary conditions)
The actual data is first fetch into frame_buffer which is a 32 bit register inside ALU (actually 16/32/64 bit depending on type of CPU) and then rest of the processing is only inside ALU's registers only.

Only when you exhaust the frame_buffer, do you really flush that bit-buffer and get another integer from memory.

This is how it gets its efficient implementation. Hope this answers your question.


Decompose video frames into slices I am working on a video motion feature implementation. I am implementing a work in a research paper. What I want to implement is this:

  The frames contained in a video clip can be thought of as a volume obtained by considering all the frames in time. This volume can be decomposed into a set of two 2-D temporal slices I(x,t)and I(y,t).

What's the methodology of implementing this i.e slicing the volume of frames in time?
Relevant link of paper: http://www.cs.cmu.edu/~yaser/RasheedSheikhShah_CVST_2005.pdf

Each frame of your movie is a 2D plane indexed by $x$ and $y$.
When you consider several consecutive frames, you can stack them in order to have a 3D cube, where each slice of the cube is a frame.
This space-time volume is now indexed by $x,y$ and $t$.

You can manipulate this volume as any 3D volume, just replace i nyour mind $t$ by $z$ if you find it confusing.

Once you have this volume, if you change the normal axis used to define the 2D slices (remember it was the $t$ axis initially), you obtain the different slices that you want:

choosing the $x$ axis as normal, you obtain 2D $(y,t)$ slices
choosing the $y$ axis as normal, you obtain 2D $(x,t)$ slices

By the way, the classical term for these kind of slices is epipolar-plane image analysis if you are loking for more references.


How do wavelet levels depend on one another? I am working with wavelets for my thesis, and I would like to know if there is an intuitive dependency between the levels of a DWT.

To make it clear: The DWT performs a convolution, and then the result is downsampled, right?

Like in this example of a convolution below.

Now when we downsample dyadically, we throw away every second value, right?
Like this:

Can we say that 2 adjacent values in the original are projected onto 1 value in the transform?

For example, the first value in the first level of the DWT would then be "responsible" for values 1 and 2 in the original.
So that we would know where e.g. a pattern in a high DWT level actually comes from.

Does this kind of projection exist?
If it does not work, why not?
Is the convolution/downsampling performed differently than I
imagined?

I may be a little confused by your question (intuitive dependency?) but I think these two sources will help you wrap your head around the relationship between levels in a DWT:

https://mil.ufl.edu/nechyba/www/eel6562/course_materials/t5.wavelets/intro_dwt.pdf
http://users.rowan.edu/~polikar/WTpart4.html

I think one way to think of it is because a DWT can be inverted, you can take a signal and perform a DWT on it and then you can reverse that DWT to get the original signal. So there's a clear dependance between levels in a DWT. I don't expect this answer to be accepted but i hope it can point you in the right direction to what you're looking for.

Small quote from link #1

"In wavelet analysis, the Discrete Wavelet Transform (DWT) decomposes a signal into a set of mutually orthogonal wavelet basis functions. These functions diﬀer from sinusoidal basis functions in that they are spatially localized – that is, nonzero over only part of the total signal length. Furthermore, wavelet functions are dilated, translated and scaled versions of a a common function φ, known as the mother wavelet. As is the case in Fourier analysis, the DWT is invertible, so that the original signal can be completely recovered from its DWT representation.
"


Impulse response of a continuous time system Question

  The impulse response of a continuous time system is $h(t) = e^{-t}u(t)$. When $x(t) = u(t - 1)$ is input to this system what is
  the value of the output signal $y(t)$ at $t = 2.34$.

My Attempt

For convolution $y(t) = h(t) \star x(t)$
\begin{align}
x(2.27) &= u(2.34-1) = 1\\
h(2.27) &= e^{-2.34}\cdot u(2.34) = 0.09632763823\\
\Longrightarrow y(t) &= 1 \cdot 0.09632 = 0.096
\end{align}

But the answer is wrong. Not too sure where my mistake is. Could someone correct my error?

Second Attempt
\begin{align}
y(t) &= ∫x(\tau)h(t-\tau)d\tau \quad\text{with limits negative and positive infinity}\\
&= ∫u(\tau-1)e^{-\tau}u(t-\tau)d\tau \\
&= ∫e^{-\tau}u(\tau-1)u(t-\tau)d\tau \\
&= ∫e^{-\tau}d\tau \quad \text{with limits 1 to 2.34}\\
&= 0.2715518029
\end{align}

Though the answer still appears to be wrong.

\begin{align}
y(2.34) &= h \star x\,\big\vert_{t=2.34}\\
&= \int_{-\infty}^\infty h(\tau)x(2.34-\tau)\,\mathrm d\tau
&\scriptstyle{\text{from the definition of the convolution integral}} \\
&= \int_0^\infty e^{-\tau}x(2.34-\tau)\,\mathrm d\tau &\scriptstyle{\text{because } h(\tau)=0\text{ for }\tau < 0}\\
&= \int_0^\infty e^{-\tau}u(1.34-\tau)\,\mathrm d\tau &\scriptstyle{\text{because }x(2.34-\tau) = u((2.34-\tau)-1)}\\
&= \int_0^{1.34} e^{-\tau}\,\mathrm d\tau
&\scriptstyle{\text{because }u(1.34-\tau) = 0 \text{ for }\tau < 1.34}\\
&= 1-e^{-1.34}
\end{align}


Real digital filter property I am a beginner to study about the filter notion and property

Being a real digital filter,
(here "real filter" I means that its impulse response is real-valued)

this formula is established.
But I have no idea how to prove it

$$ |H(\pi +w )|=|H(\pi -w)|$$

What should be $H(w)$ or $|H(w)| $ for generalized proof procedure?

HINT:

From the definition of the DTFT

$$H(\omega)=\sum_{n=-\infty}^{\infty}h[n]e^{-jn\omega}\tag{1}$$

derive the following facts:

$H(\omega)=H(\omega+2\pi)$
$H(\omega)=H^*(-\omega)$ for real-valued $h[n]$

where $*$ means complex conjugation. Combine these two results to prove the given equation.


Adequate representation of frequency domain amplitude/magnitude of FFT of a signal I'm quite new to the subject and am having fun playing around with the FFT. What I am currently doing is trying to sample an audio signal and display its frequency spectrum at the same time. This works alright with what I am doing so far, but there are still a few missing puzzle pieces.

This is what I have implemented so far:

I get the raw audio data (i.e. my signal) and sample it with a sampling frequency of 44.1 kHz.
For each time instance of the audio file where I update the spectrum (e.g. every 1/24 th second), I grab a certain amount of samples from my sampled signal around the time instance, apply a window function to it (Hann window) and then I transform it using the fftw library.
Then I look at what I get. Since all my input data is real, I get n/2 + 1 complex values, am I correct in that assumption?
Now, I'd like to display this spectrum. I calculate the magnitude of each of the n/2 + 1 vectors using the l2-norm and try to display that.

The last part is where I'm note sure of how to proceed. First, is this magnitude I calculate called the amplitude in frequency domain or am I mixing stuff up?
As far as I know, the amplitude of the original signal in the window I transformed is somehow spread across my spectrum (i.e. my frequency bins).
What I am looking for is a nice mapping in order to display these frequency bins nicely (in fact...I always group a few together to get larger frequency bands).
Here is what I am trying to do:
For my visualization, I have 12 bands and each has 15 discrete levels. So effectively what I am doing so far is trying to map the maximum of all bins in a frequency band into {0,1,...15}.

I have just been playing around so far, trying things like logarithms and linear mappings, but none seem be giving me the results I am expecting. For example, the lower frequencies might have a very high amplitude where higher ones are comparatively low, even though the corresponding audio signal would have me expecting higher magnitudes for the higher frequency bands.

So my main question, could this be mainly due to the fact that I am always picking the maximum? Would a mean value over all bins in a band be better?

My gut tells me thats just the tip of the iceberg. A linear transformation is, for example, difficult to adjust since I don't know what the maximum frequency magnitude is when I am calculating my spectrum on the fly.

I would appreciate it if some of you Gurus could help me learn what I am missing and, perhaps, tell me if there is anything horribly wrong with what I am doing up to the point where I want to visualize the calculated spectrum.

Cheers!
Brick

The correct way to group multiple bins together is to multiply each complex fft bin output by its complex conjugate (which gives the bin power) then add all the bin powers together and divide by the number of bins in the group. If you want to display in db (which is the conventional approach) then take 10*log10() of the result. Negative db values are normal and expected.

Note that if your audio signal is music, it is quite normal that the spectrum falls off with frequency.

Also from your description it sounds like you are not using a 2 ^n number of points. It's best to stick with 2^n data buffer sizes to maximize the efficiency of the fft computation.

Bob


FIR filter design for response inside unit circle I would like to design an FIR filter such that its Z-transform has a certain profile in certain regions.

For example, if I'd like to have an FIR filter that nulls the decaying exponential $\left(\frac{1}{2} \right)^n$, I can simply put a zero at $z = \frac{1}{2}$.

My question is, is there a more systematic way of designing such filter for arbitrary Z-transform response? Something similar to the Remez algorithm or maybe using convex optimization? So far all the papers I've looked at online looks at responses on the unit circle only.

Note that an FIR filter's transfer function is completely determined by its behavior on the unit circle (or any other closed contour). So if you were to define a desired $\mathcal{Z}$-transform (instead of a desired frequency response on the unit circle), you could just evaluate this desired transfer function on the unit circle and use one of the many standard design procedures.

If it were convenient for you to specify a desired response on a circle in the $z$-plane with radius $r$, you could also just use a standard frequency domain design procedure, and obtain your final coefficients by computing

$$h[n]=r^{n}\tilde{h}[n]$$

where $\tilde{h}[n]$ are the coefficients obtained from the standard design procedure.

Of course you can also directly specify the desired zeros in the complex plane, as mentioned in Jason's comment. Note that you can't specify the poles because for a causal FIR filter they must lie in the origin of the $z$-plane.


discrete wavelet transform matrix for vectorized image Yesterday I asked about how to extract 2D DFT matrix for a vectorized image. Today my question is how can I extract 2D DWT matrix for a vectorized image.
Fourier transform have this property that rows of the image are transformed first, than columns are transformed. Is there a similiar property for discrete wavelet transform such that 1D DWT matrix can be utilized for 2D DWT?

For only one level on both direction, the classical DWT schemes are applied separately (rows then columns, or colums then rows).

For more levels, there exist two schemes:

one interleaving decompositions: row, col, row, col, etc. (also called non-separable/square/non-standard/isotropic/Mallat),
one on all rows, then all colums (also called separable/rectangle/standard/anisotropic/tensor).

The second one only deserves the "fully separable" mention.

References to related questions:

2D DWT Image Issue
What is the correct order of operations for a 2D Haar wavelet decomposition?
2D DWT computation order


Can we decompose a correlated random variable into correlated plus uncorrelated? Assume a random sequence
$
Y_k  = X_k  \cdot H_k  + W_k
$
 where $X_k$ is a deterministic quantity, $H_k$ is a correlated random sequence, $W_k$ is Additive White Gaussian random sequence.

Can we decompose/rewrite the above sequence into a form like below?
$
Y_k  = X_k  \cdot \hat H_k  + X_k  \cdot \tilde H_k  + W_k
$

where $\hat H_k $ is deterministic/white and $\tilde H_k$ is correlated random - something like that?

Using Karhunen-loève expansion (principal component analysis) on a centered stochastic process $X_k$ you will get :
$$
X_k=\sum_{i=1}^{N}Z_ie_i(k) \quad k\in[1,N]
$$
Where $N$ is the process's length and $e_i$ are the $N$ eigenvectors of length $N$ from ACP. The $Z_i$ will be uncorrelated.

I'm not sure it is answering your question, I advise you to ask this on cross validated.


DFT of a complex sinusoid I'm attending this course (Coursera: Audio Signal Processing for Music Applications) in which the professor derives a general equation for Discrete Fourier Transform (DFT) for a complex sinusoid. The following is the screenshot of the slide he used:

While the derivation is totally fine, I'm having trouble understanding the concluding statement

  if $k \neq k_0, denominator \neq 0$ & $numerator = 0$
  thus $ X_1[k]=N$ for $k=k_0$ & $X_1[k]=0$ for $k \neq k_0$

I did try value substitution and all, but all that can't seem to justify the statement. My understanding is that if $k=k_0$, both the numerator and the denominator would be zero and there would be no way the result is $N$ and I have no clue about the first part of the statement either. Or I'm missing something here (or forgotten some basic school math here). What's going on here?

The last expression (sum of  a geometric series) is a common abuse of notations: it should have been:

$N$ if $k= k_0$
$f(r)=\frac{1-r^N}{1-r}$ with $r=e^{-j 2\pi(k-k_0)/N}$ if $k\neq k_0$

Indeed, as you correctly remarked, numerator and denominator would vanish for $r=1$ (or $k= k_0$), so the fraction is not "theoretically" defined. However, it is consistent to the limit, as $f(r)\to N$ as $r\to 1$  since $e^{-j 2\pi(k-k_0)/N}\to 1$ when $k \to k_0$ (allowing real $k$).

Mathematically, the fraction gives one single expression, defined by continuity at $k= k_0$. It is a bit like saying that $\frac{r^2-1}{r-1}$, not defined at $r=1$, is in some way equivalent  to $r+1$ (using $r^2-1=(r-1)(r+1)$  and wrongfully simplifying the fraction. This is not correct, but makes sense.


what does "derivative" means in image processing? In image processing and especially edge detection, when we apply sobel convolution matrix to a given image, we say that we got the first derivative of the input image, and when applying the laplacian matrix to the initial image we say that we got the second derivative.

Taking into consideration that in both cases we applied the same operation (img * matrix), then why the first operation gave us the first derivative, and the second one gave us the second derivative and not the fifth ? is the values of the mask (convolution matrix) which decide which derivative we got ? or I'm already mistaken in my description ?

Look at the numbers in the filter kernel in just 1 dimension for simplicity. For a Sobel and  Prewitt matrix you have something that roughly looks like this

[-1,0,1].

Convolving this with your image basically computes the difference between the pixel values of the neighboring pixels. You apply 0 to the current pixel, 1 to the pixel on the right and -1 to the pixel on the left. This gives a first order difference:

next pixel - previous pixel,

i.e. first derivative.

But now look at a Laplacian operator. It looks something like [1, -2, 1]. This computes the difference of differences. To see how, note that

[1,-2,1] corresponds to

next - 2 x current + previous

i.e.

next - current - current + previous

i.e.

(next-current) - (current-previous)

Now notice how this is a diference of differences. (next - current) is like a 1st derivative. (current - previous) is like 1st derivative. Their difference is like a 2nd derivative.


simple iir filter implementation in C++ I am programming the simplest way of implementing and IIR filter. However, I am struggling with the use of buffers. Do any of you have any simple implementation of an IIR filter implementing manually the convolution? I would like to have something simple, implemented on a function or something like that on C/C++.
thanks in advance.

In my opinion, you could use a couple of easiest and more optimized solutions, unless you need only that filter and you want to avoid third party or extra library/files:

SigPack: http://sigpack.sourceforge.net/index.html
Kfr: https://www.kfr.dev/

Cheers.


Time delay estimation using ultrasound I have a sound source that generates an ultrasonic pulse. I have two microphones seperated by a distance which are at a distance from the source. I want to find the time delay of arrival of the ultrasonic pulse at the mics. I am filtering out all the frequencies other than my source frequency from the outputs of the mics. Due to noise and other reasons, the shape of the pulse received at each mic is not same. I am using cross-correlation to find the time delay but the results are inconsistent. I observed that when I use low frequencies (100~300hz), results are consistent. But I want to find using ultrasonic sound. I am using the below code for finding time delay

[xc,lags]=xcorr(x1,x2);
[m,i]=max(xc);
timeDelay=lags(i);

Can anyone please tell me what is it that I am doing wrong and how to correct it? I posted a similar question here, but couldnt understand the provided answer. In some articles, it is mentioned that cross-correlation works only for cases when the signal shapes are identical. If so, how can I find the delay? I was also wondering if the voice tracking cameras used in video conferencing use some technique other than cross-correlation.

The transient edge or pulse shape of a tone burst represents a modulation.  The resulting modulated signal will contain frequency components higher than the tone frequency.  Thus you may need a sample rate much higher than just over 2X the tone frequency to reproduce a waveform accurately enough for cross-correlation to match the edges of the pulse envelopes.  This is likely true for your low tones and for voices, where the sample rate can be 10X to 200X higher than the tone or forment frequencies.

Filtering out frequencies just below the tone carrier can also remove any lower sidebands of the modulation, thus likely increasing the error in envelope matching by cross-correlation.


How $\Delta f=\frac{1}{2T}$ satifies the orthogonality condition? Let
$$s_{ml}(t)=\sqrt{\frac{2E}{T}}\exp(j2\pi\Delta fmt)$$
where $T$ is the time-period of signal
$\Delta f$ is the frequency spacing
The text says that two signals $s_{ml}(t)$ and $s_{nl}(t)$ are orthogonal to each other for $\Delta f=\frac{1}{2T}$

Orthogonality means that $$<s_{ml}(t),s_{nl}(t)>=0, \text{when}~m\neq n$$
$<>$ is inner product

On further simplifying :
$$\int_{0}^{T}\frac{2E}{T}\exp(j2\pi\Delta f(m-n)t)~\text{d}t$$
which further simplifies to
$$\frac{2E \sin(\pi T \Delta f (m-n))}{\pi T \Delta f (m-n)} \exp(j\pi\Delta f(m-n)T)$$
which is precicely
$$2E~\text{sinc}(T \Delta f (m-n))\exp(j\pi\Delta f(m-n)T)$$
As per my understanding the aforementioned term should be zero for $\Delta f=1/T$ and not $\Delta f=1/2T$.
Where I am wrong?

You're absolutely right, it must be $\Delta f = 1/T$, since this is the base frequency of the fourier transform of periodic functions with period $T$.

Otherwise, it's clear for $m=3$, for example, that the exponential would become

$g(t)=\exp\left(j \frac{2\pi}{2T}\cdot 3 \cdot t\right)$

So it does not obey periodic boundary conditions at the right boundary, i.e. at $t=k\cdot T$ because

$g(T)=\exp(3\pi j)=-1 \not= 1 = g(0)$

So the problem statement in the book is wrong, but seems to come to the right result with either wrong calculus or by already knowing the right result. Maybe the author confused $\Delta f$ with the Nyquist frequency or something. Or it was just a typo in the book.

Remark: the sinc function has zeroes for all integers $n,m$ except where $n=m$, which is kind of a pitfall of $sinc(x)$, because the limit for $x\to 0$ is one, not zero as the occurrence of the sine in the numerator alone would imply. So the orthogonality condition that you have mentioned in the question could be extended to an orthonormality condition

$$<s_{ml}(t),s_{nl}(t)>=\delta_{mn}$$

which actually also holds for $n=m$, where it presents the normalizability of the basis functions (requiring appropriate normalization factors, of course).


FIR filter design and implementation - sample rate and number of taps I have been designing and using IIR filters to process audio for a while now but want to design some FIR filters for delay equalization. I have done some preliminary work using the following approach:

I start with a desired "group delay EQ" curve. This starts at zero, and then has increasing delay at higher (audio) frequencies.
I calculate the corresponding phase lag by integrating, using the definition of group delay. This is still on the audio log-frequency grid.
In the next step, I interpolate the phase onto the equal-linear-spacing frequency grid that will be used for the discrete transform processing. I use 32768 samples for this step.
I then generate the real and imaginary components of the complex number that results for this phase angle, with the magnitude of the response equal to 1.0. I make sure that H(k) = H*(N-k).
I take the inverse DFT of the real and imaginary components to transform the response to the time domain. Although the output produces both a real and an imaginary component the imaginary component is essentially zero. The impulse is present in both short positive and negative times (at the beginning and end of the data set) but only spans about 20% of the full data set.
I set the imaginary component values all to zero. I manipulate the real component of the time domain response (the impulse), shifting it "to the right" to move the part of the impulse that is at small negative times (at the "end") to positive time. I then apply a rectangular window to the data, setting everything outside the window to zero. From what I understand, the remaining non-zero values are the FIR filter coefficients.
To find out what the filter response will look like I perform the forward DFT on the shifted, windowed impulse to bring the data back into the frequency domain. I then calculate and plot the filter magnitude (SPL), the phase, and finally the group delay.

The above procedure seems to give me the FIR filter response that I want to achieve. I can see the tradeoffs in the frequency domain when I set fewer or more taps to zero. So my approach seems to be working. Now on to the implementation and questions:

Because I am new to FIR filtering, I am confused about the data sample rate and number of coefficients calculated above and how that factors in when I want to implement the FIR filter. I may either write my own code or use some existing code or package to perform the correlation between filter and audio data stream.

Questions:

A. How will the number of bins that I use in the inverse DFT (step 5 above) influence the resulting filter? Will using more bins help to fit sharper features and vice versa? Is the number of bins usually optimized when designing a filter to keep the calculation overhead low when the filter is implemented?

B. I performed all the steps above assuming a sample rate of 48kHz. The only time that the sample rate came into play was when I calculated the frequencies of each bin (e.g., df) for the iDFT and DFT calculations. Given my final set of non-zero filter coefficients (obtained in step 6 above) is the filter simply (naively) calculated as y(n) = SUM_over_m[h(m)*x(n-m)] where 'x' is a sample in my audio data stream at the same 48kHz sample rate? That seems too easy/simple... (I understand that a practical FIR filter implementation would not carry out the calculations in this way, but would use FFT, etc.)

What really puzzles me is this: any number of coefficients can be set to zero to approximate the filter, yet you may vary the number of bins in the iDFT/DFT operations (change resolution) and then truncate the resulting filters to the same number of taps. Is it correct that the coefficients would be different even if they are the same in number, and that the coefficient values are the only thing that is "encoding" the sample rate used in the filter design and the resolution that was chosen???

Please help me understand these concepts a little more clearly. Thanks!

I'm unable to visualize the real-valued inverse-DFT FIR filter coefficients you describe in your Step 5. But no matter. My answer for A: When you use a large number of DFT bins to yield fine granularity of your frequency-domain representation of your desired frequency response, the greater will be the number of nonzero coefficients remaining after you window the output of your inverse DFT. And the larger the number of nonzero coefficients, the greater will be the number of filter computations needed per filter input sample. Only you can decide how many multiplications and additions you're willing to perform per filter input sample. By the way, I'm wondering how many nonzero coefficients you're contemplating at this time.

My answer for B: I'm not sure what the asterisk symbol means in your y(n) = SUM_over_m[h(m)*x(n-m)] expression. In any case, you typically implement your FIR filter by convolving your filter's input sequence with your filter's nonzero coefficients. If your have 100 nonzero coefficients then you must perform 100 multiplies and 99 additions per filter input sample. (Depending on your computing hardware, there may be a way to reduce the number of multiplies per input sample to 50. Ask me how.) Now if you have many hundreds of nonzero coefficients then in that case, yes, you'll use what's called the "fast convolution" filtering method which involves both forward- and inverse-DFTs. But CA Charlie, all of this discussion depends on whether your filtering is "real-time" of "offline." As for the 48 kHz data rate, it affects two things: (1) it means your input signal's spectrum can contain no spectral components greater than 48/2 = 24 kHz; and (2) in a real-time implementation you must perform all your "per sample" filtering computations in less than 1/48,000th of a second. (Sorry CA Charlie, I don't understand the question in your final paragraph.)


How to transform SoundWave Picture to sound? I would like to recreate a sound from a picture of a soundwave.

But i see a lot of theorical question.

And im not sure its possible....

I mean, it is possible to create a sound...

But to make the sound as close as possible of the original one, sound a bit harder.

Most of the time when you see soundwave as art...

It's really time compressed.

I would love to have a tatoo of the soundwave of someone important to me that say a word to me.

And then with my smartphone point the camera at the tatoo and hear it again.

But i have a lot of question about this.

First how can i get the periodicity?
Per pixel i will probably got a lot of periode.
How can i decide the lenght of my whole song?
I have a sound wave but no legend. so it can be 1s or 160s compressed i can never know....

Let's say i cheat and decided my application always says from start to end of my soundwave i have 1s (1s is enough for a 'I love you' for exemple)

Ok so i have my length total so with open cv i can check my image and reconise bar height every pixel.
And then with my pixel length and the fact i know it's one second i can start to generate sound.

But how can i decided.... the corresponding value to the Heigh in pixel? I don't have frequency value correspondance. 12px heigh can be 20Hz as 20 000Hz

Ok let's say i define a value to the max height of my sound wave to the real value of my sound.
Cause i use a sound i choose.
So i can set a reference value and use it to get other value.

I now parse a picture to get a compressed sound and i have some value of Hz in an intervall of second.

Lets say i have a 10cm tatoo of a 1s sound.
each cm is 0.1s

And in each cm i can make like... 10 bar maybe 20 (one bar is 1mm or 0.5mm possible in tatoo right?)
Let say 20 ...

so i have info of Hz every
0.005s

and let say i decide i start up and and bottom
i hav two info in 0.005s the up and bottom one
So i have one info for 0.0025s

That make a sound at 400Hz...
And... classical sound are know to be around 44 100Hz....

Will the sound i produce be at little close to the real one?
Or does it will only be a little noise?

(English is not my mother tongue feel free to improve i will try my best to reformule if you didn't understand something)

I don't think "I love you" is doable with a sound wave image, because that's already very roughly about 200 cycles of the fundamental frequency of the voice, and you need also the harmonics so that it sounds right and can be understood. Too much horizontal resolution is needed to pack that into a tattoo.

You could probably pull it off as a spectrogram.

I've seen ads for a commercial service that presumably stores at a server a database of sound clips each associated with a sound wave image. Then their phone app sends the tattoo photo to the server and plays back the sound clip associated with the image that most resembles the tattoo photo. Technically that's not far different from having the sound clip number tattooed as a bar code.


Implement Fast Hartley Transform I'm trying to implement a simple FHT on .NET platform myself. I follow this document: Optimized fast Hartley transform for the MC68000 with applications in image processing.
In the page 20th, equation (20):

  H(k) = He(k) + [Ho(k)cos(2πk/N) + Ho(N/2-k)sin(2πk/N)] H(k+N/2) = He(k) -
  [Ho(k)cos(2πk/N) + Ho(N/2-k)sin(2πk/N)]

  He(k) is the N/2 point DHT of the even indexed elements of
  H(k) and
  Ho(k) is the N/2 point DHT of the odd indexed
  elements of H(k)

$k$ is zero index, so where is the item $N/2$ if $N/2$ is the length? I know I must have misunderstood something, but I realy don't know what exactly the Ho and He are. I have succeeded implement the FFT equation (in the page 17th). But the FHT, I can not.
Can somebody please help me explain the above equation details?

First of all, I'm pretty sure that $H_e(k)$ is the $N/2$ point DHT of the even elements of the original sequence to be transformed, not of $H(k)$. The same is true for $H_o(k)$. As for the index $k$, when used with $H_e(k)$ and $H_o(k)$, it is evaluated modulo $N/2$, so it ranges from $0$ to $N/2-1$. So the index $k=N/2$ is mapped to the index $0$ when evaluating $H_e(k)$ and $H_o(k)$.


Melody/Vocal extraction I'm new to the whole DSP area, but I've built a program which plays music and does FFT to produce some cool visuals for the music in realtime.

The FFT algorithm I chose was from a C# library called NAudio. I applied a hamming window before shipping it off to the FFT, from which I get an array of complex. Nothing special. Works fine.

But now I want to display this data to the user, so I converted the bins to magintudes following this formula: magnitudes[i] = sqr(real[i]² + imaginary[i]²)

This with some amplification yields a pretty good result to display for the user, but I feel that most of the levels being displayed are just random noise, I want to show the user what he/she "hears", where the melody or vocal are being amplified OVER the other background music.

The thing I'm unsure of is how to approach this problem. What can I do to filter out everything else and only have the melody/vocals left. I can understand there's not 100% clean way you can do this, but what method could give me a decent result?

Thank you

I remember reading a similar question (either here or on stack overflow) abouylt melody extraction. In any event someone did create such a project as a PhD thesis (luckily I bookmarked the site). He even offers his work (for research/personal/non-commercial use) for windows mac and linux http://www.justinsalamon.com/melody-extraction.html#demo if you could somehow integrate this to your application you can use the extracted melody to display the visuals, but still just play the original song from speakers/headphone


Recursive filter gain I implemented a filter similar to the one below, within a FPGA. Two substractors and an accumulator.

The recursive equation is given by:

My inputs are exponentially decaying pulses of the form  sampled every n

I want to express y[n] as a function of the pulse amplitude A (which I need to retrieve), and the constants , L, and K.

Any hints would be greatly appreciated

Edit:

I followed Andy Walls' advice and this is what I get

$$H(z) = \frac{(1-Z^{-K})(1-Z^{-L})}{1-Z^{-1}}$$
$$X(z) = \frac{A}{1-aZ^{-1}}   $$
$$with: a = e^{\frac{-nT}{\tau}}$$
$$Y(z) = H(z).X(z) = \frac{A}{1-aZ^{-1}} .  \frac{(1-Z^{-K})(1-Z^{-L})}{1-Z^{-1}} $$
Using partial fraction expansion I try to find B and C so that
$$ Y(z) = A.(\frac{B}{1-aZ^{-1}} + \frac{C}{1-Z^{-1}})$$
I get $$ B = \frac{A.(1-a^{-K})(1-a^{-L})}{(1-a^{-1})}; C=0$$
So $$Y(z) = \frac{A.(1-a^{-K})(1-a^{-L})}{(1-a^{-1})(1-aZ^{-1})}$$

I don't know if I am right at this point, but this is where I get stuck. This filter is trapezoidal, but because the input is an exponential decay, I would expect at least 2 exponential terms corresponding to the trapeze rise and fall times.

Am I on the right path?

Take the Z Transform of your difference equation and perform algebraic manipulations to find $H(z)=\dfrac{Y(z)}{X(z)}$ as a rational function of $z$.

Take the Z Transform of your input $x(nT) = x[n]$ to find $X(z)$.

Multiply to get $H(z)X(z) = Y(z)$.

Find the inverse Z Transform of $Y(z)$ to find $y[n]$.  (This will be the difficult part, as you may need to use partial fraction expansion and other mathematical techniques so you can look up inverse transform in a table.)

These tables will help:
http://lpsa.swarthmore.edu/LaplaceZTable/LaplaceZFuncTable.html

Update in response to comment

Here's what I get.

Finding the transfer function:

$$\mathcal{Z}\left\{y[n]\right\} = \mathcal{Z}\left\{y[n-1] + x[n] -x[n-K]-x[n-L]+x[n-K-L]\right\} $$

$$Y(z) = Y(z)z^{-1} + X(z) - X(z)z^{-K} -X(z)z^{-L} + X(z)z^{-K-L}$$

$$\begin{align*}H(z) = \dfrac{Y(z)}{X(z)} &= \dfrac{1-z^{-K}-z^{-L}+z^{-K-L}}{1-z^{-1}}\\
\\
&= \dfrac{\left(1-z^{-K}\right)\left(1-z^{-L}\right)}{1-z^{-1}}\\
\\
&= \left(1 + z^{-1} + \dots + z^{-K+1}\right)\left(1-z^{-L}\right)\\
\\
&= 1 + z^{-1} + \dots + z^{-K+1} -z^{-L}-z^{-L-1}- \dots -z^{-K-L+1}\\
\\
\end{align*}$$

Transforming the input:

$$X(z) = \mathcal{Z}\left\{x(nT)\right\} = \mathcal{Z}\left\{Ae^{-\frac{nT}{\tau}}u[n]\right\} = A\dfrac{z}{z-e^{-\frac{T}{\tau}}}= A\dfrac{1}{1-e^{-\frac{T}{\tau}}z^{-1}}$$

Filtering in the Z domain:
$$\begin{align*}Y(z) &= H(z)X(z)\\
\\
&=\dfrac{1-z^{-K}-z^{-L}+z^{-K-L}}{1-z^{-1}} \cdot A\dfrac{1}{1-e^{-\frac{T}{\tau}}z^{-1}}\\
\\
&=A\left[\dfrac{1-z^{-K}-z^{-L}+z^{-K-L}}{\left(1-z^{-1}\right)\left(1-e^{-\frac{T}{\tau}}z^{-1}\right)}\right]\\
\\
&=A\left[\dfrac{C}{1-z^{-1}}+\dfrac{B}{1-e^{-\frac{T}{\tau}}z^{-1}}\right]\\
\\
\end{align*}$$

Performing partial fraction expansion:

$$\begin{align*}1-z^{-K}-z^{-L}+z^{-K-L} &= C\left(1-e^{-\frac{T}{\tau}}z^{-1}\right) + B\left(1-z^{-1}\right)\\
\\
\left(1 + z^{-1} + \dots + z^{-K+1}\right)\left(1-z^{-L}\right)\left(1-z^{-1}\right)&=C\left(1-e^{-\frac{T}{\tau}}z^{-1}\right) + B\left(1-z^{-1}\right)\\
\end{align*}$$

so
$$\begin{align*}C &= 0 \\
\\
B &= 1 + z^{-1} + \dots + z^{-K+1} -z^{-L}-z^{-L-1}- \dots -z^{-K-L+1}\\
\end{align*}$$

Continuing filtering in the Z domain:
$$\begin{align*}Y(z) &=A\dfrac{1 + z^{-1} + \dots + z^{-K+1} -z^{-L}-z^{-L-1}- \dots -z^{-K-L+1}}{1-e^{-\frac{T}{\tau}}z^{-1}}\\
\\
\end{align*}$$

Since there is no information on the values of $K$ and $L$, there is not really any more simplification we can do in the Z domain.

Taking the inverse Z-Tranform:

$$\begin{align*}y[n] &= A\left[e^{-\frac{nT}{\tau}}u[n]+e^{-\frac{(n-1)T}{\tau}}u[n-1]+\dots+e^{-\frac{(n-K+1)T}{\tau}}u[n-K+1]-e^{-\frac{(n-L)T}{\tau}}u[n-L]-e^{-\frac{(n-L-1)T}{\tau}}u[n-L-1]-\dots-e^{-\frac{(n-K-L+1)T}{\tau}}u[n-K-L+1] \right]\\
\\
&= A \sum_{k = 0}^{K-1} e^{-\frac{(n-k)T}{\tau}}u[n-k] - e^{-\frac{(n-k-L)T}{\tau}}u[n-k-L]\\
\\
\end{align*}$$


High frequencies disappear when applying discrete wavelet transform Trying to decompose and reconstruct a signal using a to some extent self-made implementation of DWT for some reason fails. The result looks highpass filtered and/or shifted. I wanted to write the code in Python not using the PyWavelets package (except for getting the wavelet object).

import numpy as np
from scipy import interpolate, signal
import pywt
import matplotlib.pyplot as plt

# Test signal and wavelet
sig = np.zeros((10*512))
for i in range(256):
    sig += i * np.cos(2*np.pi*i*np.linspace(0, 10, 10*512)+i**2*2*np.pi/((i+1)/1024))
w = pywt.DiscreteContinuousWavelet('db6')

I convolved the test signal with the decomposition filters from the wavelet object and used scipy's signal.decimate(...). As far as I understood, the decomposition filters do not lowpass filter the signal (or do they and is this where I went wrong?). However, a simple downsampling instead of decimation did not yield better results.

# Downsampling by decimation, decomposition
coefs_own = []
cA = np.copy(sig)
for i in range(1):
    cA_last = cA
    cA = np.convolve(cA_last, w.dec_lo, 'same')
    cD = np.convolve(cA_last, w.dec_hi, 'same')
    cA = signal.decimate(cA, 2)
    cD = signal.decimate(cD, 2)
    coefs_own.insert(0, cD)
coefs_own.insert(0, cA)

And this is my reconstruction using linear interpolation and a convolution with the recomposition filters provided by the wavelet object.

# Helper function
def upsample(data):
    x = np.array(range(data.shape[0]))
    x_double = np.linspace(x[0], x[-1], x.shape[0]*2)
    f = interpolate.interp1d(x, data, kind='linear')
    data_upsampled = f(x_double)
    return data_upsampled

# Upsampling by interpolation, reconstruction
cA = coefs_own[0]
for i in range(1):
    cA_up = upsample(cA)
    cD_up = upsample(coefs_own[i+1])
    cA = np.convolve(cA_up, w.rec_lo, 'same') + np.convolve(cD_up, w.rec_hi, 'same')
sig_rec_own = cA

For comparison, see the multilevel decomposition by PyWavelets, which does a perfect reconstruction.

# PyWavelets
coefs = pywt.wavedec(sig, w, mode='zero', level=1)
sig_rec_pywt = pywt.waverec(coefs, w, mode='zero')

# Plot both decomposition versions
fig, axes = plt.subplots(2, 1)
axes[0].plot(coefs_own[0])
axes[0].plot(coefs_own[1])
axes[0].set_title('Own')
axes[0].legend(['Approx.', 'Detail'])
axes[1].plot(coefs[0])
axes[1].plot(coefs[1])
axes[1].set_title('PyWavelets')
axes[1].legend(['Approx.', 'Detail'])

# Plot both reconstructions
fig, axes = plt.subplots(3, 1)
axes[0].plot(sig)
axes[0].set_title('Original')
axes[1].plot(sig_rec_own)
axes[1].set_title('Self-made reconstruction')
axes[2].plot(sig_rec_pywt)
axes[2].set_title('PyWavelets reconstruction')

plt.show()

Does anybody have an idea, why the decomposition and the reconstruction both obviously are erroneous? What did I miss in my implementation?

The wavelet transform is usually constructed in such a way that in the reconstruction the input signals are either convoluted with each of the downsampled filters and then summed up, or equivalently, are upsampled filling the gaps with zeros.

Upsampling with interpolation is equal to applying an additional filter $[0.5,1,0.5]$ which is a low-pass, smoothening filter. Esp. on the high-frequency or detail component this will lead to large errors.


Similarity or Relation between Walsh Hadamard Transform and Slant Transform Is there any relation between Walsh Hadamard Transform and Slant Transform? Or is there any common property or similarity?

On a 2-sample vector, or $2\times2$ pixel images, Haar, Walsh-Hadamard, or Slant are exactly the same. They are all separable, orthogonal transformations with fast $N\log N$ algorithms. According to Slant Transform Image Coding, 1974, it is superior in appearance to the above, only inferior to Karhunen-Loève,  which has no fast enough algorithms.


Nyquist rate determination and oversampling I'm new to this. I'm using ShannonGui to generate a sine wave with a $1\textrm{ Hz}$ signal frequency.

From what I can tell to determine the Nyquist rate I need to double the highest frequency, since it's a constant $1\textrm{ Hz}$, I just double it and get a Nyquist rate of $2\textrm{ Hz}$ correct?

Now I also read that the Nyquist rate is the lowest frequency that you can use as a sampling frequency to avoid aliasing, but when I use $2\textrm{ Hz}$ as the sample frequency the reconstructed wave is a straight line at 0. So does this mean the Nyquist rate is wrong or my understanding of it is wrong?

you actually have to sample at more than twice the highest frequency.  if $B$ is the bandwidth or highest frequency and $f_\text{s}$ is the sampling frequency, then to satisfy the Sampling Theorem (for normal baseband sampling):

$$ 2B < f_\text{s} $$

it is not

$$ 2B \le f_\text{s} $$

sampling at exactly twice the frequency of a sinusoid is "critically sampling" and loses amplitude or phase.  it can't get both.  in your case, it is sampling the sinusoid when it crosses zero.


Phase vocoder, pitch shift in MATLAB and C I'm trying  to make a Phase Vocoder in MATLAB. But I'm stuck with the calculations on changing the pitch. I browsed loads of website and I simply do not understand what to do.

Currently I have in my MATLAB code a working flow where input = output

$${\xrightarrow{\rm Input}}\boxed{\text{sampling}}{\rightarrow}\boxed{\text{Hanning window}}{\rightarrow}\boxed{\text{STFT}}{\rightarrow} \cdots {\rightarrow}\boxed{\text{ISTFT}}{\rightarrow}\boxed{\text{Add}}{\xrightarrow{\rm Output}}$$

What do I have to do in the $"\cdots"$ to get the pitch shift?
With the STFT I get the FFT of 1024 sample windows with 512 overlap.

alright, conceptually, for one frame (which is centered at $n=0$ and has non-zero width of $N-1$ samples):

$$\begin{align}
 X[k] &= \mathcal{DFT}\Big\{ x[n] w(n) \Big\} \\
\\
  &= \sum\limits_{n=-\tfrac{N}{2}}^{\tfrac{N}{2}-1} \big( x[n] w(n) \big) \, e^{-j \frac{2 \pi}{N} nk} \qquad |k| \le \tfrac{N}2 \\
\end{align}$$

here, our frame hop is $\tfrac{N}2$ samples, which is half of the window width.

$w(t)$ is a continuous-time complementary window with non-zero width of $N$ samples and centered at $t=0$, like a Hann window:

$$ w(t) = \begin{cases} \tfrac12 + \tfrac12 \cos\left(\tfrac{2\pi}{N}t \right) \qquad & |t| < \tfrac{N}2 \\
 0 \qquad & |t| \ge \tfrac{N}2 \\
\end{cases} $$

the output spectrum is a scaled copy of the input spectrum.  this is what shifts the pitch.

$$ Y[k] = \begin{cases}
 X(\alpha^{-1} k) \qquad & |k| \le \tfrac{N \alpha}{2} \\
 0 \qquad & |k| > \tfrac{N \alpha}{2} \\
\end{cases} $$

where $\alpha$ is the output-to-input frequency ratio.  $X(f)$ is a continuous-frequency function that is interpolated from the discrete-frequency spectrum $X[k]$ in such a way that equality exists when the argument is an integer.

$$ X(f) \bigg|_{f=k} = X[k] \qquad k \in \mathbb{Z} $$

the discrete-time-domain output for this frame (before adjusting for stretching or scrunching the window) is:

$$\begin{align}
 y[n] &= \mathcal{iDFT}\Big\{ Y[k]  e^{j \phi[k]} \Big\} \\
\\
  &= \tfrac{1}{N} \sum\limits_{k=-\tfrac{N}{2}}^{\tfrac{N}{2}-1} \big( Y[k]  e^{j \phi[k]} \big) \, e^{j \frac{2 \pi}{N} nk}  \qquad |n| \le \tfrac{N}2 \\
\end{align}$$

and the frame of output (centered at time $n=0$) is

$$ y[n] \frac{w(n)}{w(\alpha \, n)}  $$

the $\frac{w(n)}{w(\alpha \, n)}$ factor is undoing the stretched or scrunched input window and re-applying the original window at the same scale (since we're not time-scaling and the frame hop length remains unchanged from input to output).

$\phi[k]$ is a phase adjustment (with odd symmetry in $k$) that is constant for each sinusoidal component (a la Miller Puckette) rather than a changing phase adjustment for each FFT bin (a la Portnoff).  that phase adjustment for each frequency component is what is required to make each sinusoid component continuous because of the frequency shift.  that's how the phase vocoder does glitch-free pitch shifting, even if the input is not periodic or a single harmonic note.

this is why this phase adjustment is needed.  consider a single sinusoid with normalized angular frequency $\omega_0$.

$$\begin{align}
 x[n] &= \cos(\omega_0 n + \theta_0) \\
  &= \tfrac12\big(e^{j (\omega_0 n + \theta_0)} + e^{-j (\omega_0 n + \theta_0)}  \big) \\
  &= \tfrac12\big(e^{j\theta_0}e^{j \omega_0 n} + e^{-j\theta_0}e^{-j \omega_0 n} \big) \\
\end{align}$$

now let's consider only the positive frequency complex component.

$$ \hat{x}[n] = e^{j\theta_0}e^{j \omega_0 n} $$

this has an instantaneous angle of $\theta_0$ at time $n=0$.  now the instantaneous value of this sinusoid in the center of the previous frame is

$$ \hat{x}[-\tfrac{N}{2}] = e^{j\theta_0}e^{-j \omega_0 \frac{N}{2}} $$

and the instantaneous angle is $\theta_0 - \tfrac{N}{2} \omega_0$.  and it shouldn't surprise us that halfway between the two frame centers (which is exactly in the middle of the crossfade from the previous frame to the current frame) the instantaneous angle $\theta_0 - \tfrac{N}{4} \omega_0$

now if the output spectrum is mapped from the interpolated input spectrum as

$$ Y[k] = X(\alpha^{-1} k) $$

then this component frequency at the input, $\omega_0$ gets mapped to $\alpha \omega_0$.  and the current frame output (before the phase adjustment) is

$$\begin{align}
 y[n] &= \cos(\alpha\omega_0 n + \theta_0) \\
  &= \tfrac12\big(e^{j (\alpha\omega_0 n + \theta_0)} + e^{-j (\alpha\omega_0 n + \theta_0)}  \big) \\
  &= \tfrac12\big(e^{j\theta_0}e^{j \alpha\omega_0 n} + e^{-j\theta_0}e^{-j \alpha\omega_0 n} \big) \\
\end{align}$$

and the positive frequency component is

$$ \hat{y}[n] = e^{j\theta_0}e^{j \alpha \omega_0 n} $$

and the value at the left edge of the current frame where it is joined to the right edge of the previous frame

$$ \hat{y}[-\tfrac{N}{4}] = e^{j\theta_0}e^{-j \alpha \omega_0 \frac{N}{4}} $$

now the instantaneous angle of the output sinusoid at the center of the previous frame is the same instantaneous angle of the input sinuosoid at the center of the previous frame which is $\theta_0 - \tfrac{N}{2} \omega_0$.  this makes the output sinuosoid (the positive frequency component) of the previous frame, when "justified" or expressed in terms of the time index of the current frame

$$ \hat{y}[n] = e^{j(\theta_0 - \frac{N}{2}\omega_0)}e^{j \alpha \omega_0 (n + \frac{N}{2})} $$

(note that when $n=-\tfrac{N}{2}$, which is the center of the previous frame, the instantaneous angle is $\theta_0 - \tfrac{N}{2} \omega_0$.)  now at the output sample midway between the previous frame and the current frame (which is $n=-\tfrac{N}{4}$), the sinusoid of the previous frame is at angle

$$\theta_0 - \tfrac{N}{2}\omega_0 + \alpha \omega_0 (-\tfrac{N}{4} + \tfrac{N}{2})
 = \theta_0 + (\alpha - 2)\omega_0\tfrac{N}{4}$$

but, at that same sample, the angle of the sinusoid of the current frame is

$$ \theta_0 - \alpha \omega_0 \tfrac{N}{4} $$

note that when there is no pitch shift and $\alpha = 1$, then the two angles are the same and the splice is seamless.  but when $\alpha \ne 1$, the phase of one of the frames must be adjusted to make the splice seamless.  the previous frame is already a "done deal" and will not be modified, so it's the phase of the sinusoid of the current from that must be adjusted to make the splice seamless.  that phase is modified for all DFT bins associated with this particular sinusoid as

$$\begin{align}
 \theta_0 - \alpha \omega_0 \tfrac{N}{4} + \phi[k] &= \theta_0 + (\alpha - 2)\omega_0\tfrac{N}{4} \\
 -\alpha \omega_0 \tfrac{N}{4} + \phi[k] &= (\alpha - 2)\omega_0\tfrac{N}{4} \\
  \phi[k] &= (\alpha - 1)\omega_0\tfrac{N}{2} \\
\end{align}$$

so the phase adjustment for the current frame necessary to align it with the previous frame is the normalized angular frequency of the sinusoid, $\omega_0$, times the frame hop displacement $\tfrac{N}{2}$ times the factor $(\alpha - 1)$.  you should apply same phase adjustment to all DFT bins of $Y[k]$ that are associated with this particular sinusoidal component at frequency $\omega_0$.  the values of $k$ will be around $\tfrac{\omega_0}{2 \pi} N$.

so Miller Puckette says

$$\phi[k] = (\alpha - 1) \omega_0 \tfrac{N}{2} \qquad \text{for } k \approx \tfrac{\omega_0}{2 \pi} N$$

and Portnoff would say that the adjustment would be

$$\phi[k] = (\alpha - 1) \tfrac{2 \pi k}{N} \tfrac{N}{2} = (\alpha - 1) \pi k$$

but, in my opinion, Puckette is correct and Portnoff is wrong.


Why does MATLAB's dualtree3 function not return the lowpass subband of the imaginary tree? MATLAB has a function named dualtree3 which computes 3-D dual-tree complex wavelet transform. This function only retains the last low-pass coefficients/subband of the real tree. And, it does not save the low-pass subband of the imaginary tree. How can you justify this?

Prof. Nick Kingsbury kindly provide an answer to my question!.

In 1-D, the lowpass basis functions (scaling functions) from the two
trees (a) and (b) of the dual tree WT tend to look very similar to
each other, apart from a shift of half the output sample period
between them. Hence it usually makes most sense to regard the two sets
of lowpass samples as interleaved samples from a single 2x
over-sampled real basis function.
On the other hand, the highpass basis functions (wavelets) from (a)
and (b) look approximately like a Hilbert pair of functions. so they
are most logically interpreted as the real and imaginary parts of an
analytic basis function.
Hence the over-sampling by 2 in the lowpass channel replaces the
spitting into real + imaginary channels for the bandpass channels of a
multi-scale wavelet transform.  Similar things apply when 4 trees are
used for the 2-D dual-tree WT, or when 8 trees are used for the 3-D
version.
I hope this helps.  Do refer to our 2005 paper in IEEE Signal
Processing Magazine for more details.

[UPDATE at 9/18/2020]
Since I quoted Prof. Kingsubury's answer above, I should also note one more thing. He mentioned in his email that he is not familiar with dualtree3 function in MATLAB. However, based on my question, he gave the above general answer.


Spectral density of modulated noise process Suppose we have a noise process $F(t)$.
Suppose $F$ is stationary and we know its autocorrelation function
$$\rho_F(\tau) \equiv \langle F(0)F(\tau) \rangle$$
where the average is over realizations of the process.
We then construct a new noise process $G(t) \equiv \cos(\Omega t)F(t)$.
Obviously, $G$ is not stationary.
The autocorrelation function of $G$ is
\begin{align}
\langle G(t_1) G(t_2)\rangle
&= \langle F(t_1) \cos(\Omega t_1)F(t_2) \cos(\Omega t_2) \rangle \\
&= \langle F(t_1) F(t_2) \rangle \cos(\Omega t_1) \cos(\Omega t_2) \\
&= \rho_F(t_2 - t_1) \cos(\Omega t_1) \cos(\Omega t_2) \, .\\
\end{align}

The spectral density of $F$ can be usefully defined as
$$S_F(\omega) = \int \frac{d\omega}{2\pi} \rho_F(\tau) e^{-i \omega \tau} \, . $$
However, since $\langle G(t_1) G(t_2) \rangle$ depends on $t_1$ and $t_2$ and not just their difference, we cannot easily define $S_G(\omega)$.
This is distressing because for deterministic signals we expect modulation by a sinusoid to just shift the spectral density, so we expect there should be a simple formula for $S_G(\omega)$.

What, if any, is a meaningful way to construct the spectral density of $G$?

This is a standard problem in (among others) the theory of analog or digital communication systems. If $F(t)$ is real-valued, then the standard gimmick is indeed - as already pointed out in a comment by Dilip Sarwate - to add a random phase to the carrier phase. This is not only a trick but it actually reflects reality in the sense that we often do not know (and do not care to know) the exact phase of the carrier with respect to the signal $F(t)$. So what we actually want is NOT the power spectrum of $F(t)\cos(\omega_0 t)$ (which does not exist), but the power spectrum of the process $F(t)\cos(\omega_0 t+\theta)$ with some unknown phase $\theta$, the obvious model of which is a random variable with a uniform distribution in $[0,\pi)$. This random phase is of course constant over time, but its value is chosen at random reflecting our uncertainty about the origin of the time axis. With this random phase the resulting process is wide-sense stationary (WSS) (if $F(t)$ is), and the power spectrum looks exactly as one would expect.

Now comes the interesting part. If $F(t)$ is a WSS complex-valued process and if we're interested in the power spectrum of the process

$$G(t)=\text{Re}\{F(t)e^{j\omega_0 t}\}\tag{1}$$

then we could also use the random phase trick, but under some reasonable assumptions we don't even need it. The process $G(t)$ in (1) occurs e.g. in systems using quadrature modulation. It is straightforward to show that the autocorrelation of $G(t)$ is

$$\begin{align}R_G(t,\tau)&=E\{G(t+\tau)G^*(t)\}=\\&=\frac12\text{Re}\{R_F(\tau)e^{j\omega_0\tau}\}+\frac12\text{Re}\{R_{FF^*}(\tau)e^{j\omega_0(\tau+2t)}\}\tag{2}\end{align}$$

which becomes independent of $t$ if $R_{FF^*}(\tau)=0$, i.e. if the cross-correlation of $F(t)$ and $F^*(t)$, vanishes. If $F(t)=A(t)+jB(t)$, it can be shown that this is the case if

$$R_A(\tau)=R_B(\tau)\tag{3}$$

and

$$R_{AB}(\tau)=-R_{AB}(-\tau)\tag{4}$$

So if conditions (3) and (4) are satisfied, we don't even need a random carrier phase in order to define a power spectrum of $G(t)$ as defined in (1). Conditions (3) and (4) mean that the real and imaginary parts of $F(t)$ have the same autocorrelation, i.e. the same power spectrum, and that their cross-correlation is an odd function, i.e. their cross-power spectrum is purely imaginary. This last condition also implies $R_{AB}(0)=0$, i.e the real and imaginary part are uncorrelated when sampled at the same time. It turns out that for many practical problems these conditions are at least approximately satisfied.

Given that conditions (3) and (4) are satisfied, the autocorrelation function of $G(t)$ is given by

$$R_G(\tau)=\frac12\text{Re}\{R_F(\tau)e^{j\omega_0\tau}\}\tag{5}$$

and its power spectral density is given by the Fourier transform of (5):

$$S_G(\omega)=\frac14\left[S_F(\omega-\omega_0)+S^*_F(-\omega-\omega_0)\right]\tag{6}$$

where $S_F(\omega)$ is the power spectral density of $F(t)$.


Is it valid to increase amplitude (and presumably FFT quality) by simply scaling the data? I'm using a version of "KISS FFT" by Mark Borgerding.  It accepts an array of 16-bit fixed-point input values and produces a 32-bit float result array.

I've discovered that if the input amplitudes are low many of the float result values come out zero, but if I simply scale the inputs (by, say, the factor 16) then fewer output values are zero and therefore the output seems to contain more detail.  (Not that it matters much for my purposes, but for consistency I then divide the resulting float values by the same scaling factor.)

Anyway, this seems to work, in terms of producing a result when previously I'd have just gotten a buffer of virtually all zeros, but I'm wondering if there's some reason it might not be a valid approach.

(Note that this approach means that there's a lot more "coarsness"/granularity to the data, and, in particular, the low-level noise that would normally be present is not.  I'm almost wondering if it would be wise to inject some low-level noise to replace the zero values in the input.)

This can be a valid approach. You're observing a very practical issue that arises often when using fixed-point (i.e. integer) arithmetic (although it can happen in floating-point also). When the numeric format that you're using to perform calculations doesn't have enough precision to express the full range of values that can arise from your calculations, some form of rounding is required (e.g. truncation, round-to-nearest, and so on). This is often modeled as an additive quantization error to your signal.

However, for some combinations of algorithm and rounding scheme, when the magnitude of the input signal is very low, it's possible to get what you observed: a large number of zero outputs. Basically, somewhere in the sequence of operations, the intermediate results are becoming small enough to not break the threshold required to quantize to a nonzero level. The value is then rounded to zero, which can often propagate forward to the output. The result is, as you noted, an algorithm that generates a lot of output zeros.

So can you get around this by scaling the data up? Sometimes (there are very few techniques that work all the time!). If your input signal is bounded in magnitude to a value below full-scale of the numeric format (16-bit signed integers run from -32768 to +32767), then you could scale the input signal up to more fully use the range available to it. This can help to mitigate the effects of roundoff error, as the magnitude of any roundoff error becomes smaller in comparison to the signal of interest. So, in the case where all of your outputs are getting rounded to zeros internally to the algorithm, this may help.

When can such a technique hurt you? Depending upon the structure of the algorithm's calculations, scaling the input signal up might expose you to numeric overflows. Also, if the signal contains background noise or interference that is larger in magnitude than the algorithm's roundoff error, then the quality of what you get at the output is going to be typically limited by the environment, not error introduced in the computation.


Is it normal that a DAC outputs spikes? I'm using a Raspberry Pi 2 and an AD/DA Board, with no extra circuit. Reading on the internet though, I found that DACs need a low pass filter in their output, and I don't understand if the behaviour of my DAC is normal, and I wonder if I need the low pass filter.
Its output goes down to 0V at every sample, is it normal? Is this the reason why low pass filters are used?

Thanks in advance,

Alberto

I tried with a better oscilloscope to analyze the signal, and it's all fine. I'm so sorry if I created a useless thread, but I still have learnt new terms I didn't know before.

The TI DAC8532 uses a Kelvin resistor-string (string of resistors that divide the reference voltage, with individual switches from the nodes to the output), so it is guaranteed to be monotonic. What you see does not sound normal, you must be writing zeros to it over the serial interface somehow. How did you implement the serial protocol?

The low-pass filter you put on the output is the reconstruction filter. It is used to reduce imaged frequency content (repeated frequency content above the Nyquist frequency). It's not always needed. The DAC has a built-in zero-order hold element that attenuates the images somewhat. If you are driving a system with low-pass frequency response, like a DC motor with a heavy mass attached, then you don't need the reconstruction filter.


What is the purpose of dividing an audio signal into segments and analysis each segment? I read bunch of materials for extracting feature from audio signal and they all tell me to break signal into segments, why don't we analyze all the audio signal? I don't know what are the advantages of doing that and how wide a segment should be? I only see 256 samples per frame or 512 samples per frame... what about 1028 per frame?

Analyzing signals per segments, with proper windowing, is a way to cope with non-stationary in audio samples.  With full-size analysis, features can get mixed. Segment-splitting is thus at play in many algorithms (mp3, shazam).
The length of window is often a matter of trade-offs, between data information and computing advantages:

signal sampling (window length is quite meaningless without sampling rate), with respect to the following:
analyzing or extracting informational content from the signal: various ranges of stationarity may exist in the data, or generally useful processing features,
easiness in computing: the power-of-two length you mention can be beneficial (faster algorithms like in the FFT), parallel computing, dedicated hardware, closer to real-time analysis.


Example of an LTI system with complex impulse response In general, I know that the impulse response $h(t)$ of an LTI system can be complex. However, all of physically realizable, useful systems I've come across have purely real impulse responses. I did a web search in an attempt to find a useful, non-pathological system with a complex impulse response, but was not immediately successful.

Could anyone provide an example of a system with a complex impulse response, and state what applications it has? Links to such examples would also be appreciated.

Complex band pass filters are used to get a band and its quadrature in one efficiently computed step. The most simple such design is a single-pole resonant bandpass with the discrete transfer function

$$H(z)=\frac{(1-r)\exp(i \phi)}{z-r \exp(i \phi)}$$

where $r<1$ determines the bandwidth and $\phi \in [0,2\pi[$ the band center frequency. If bandwidth and center frequency are chosen to avoid significant negative frequency amplitudes, the resulting signal is practically analytic.

Strongly bandlimited and analytic signals have a complex magnitude envelope that is well behaved and can be interpreted as the instantaneous amplitude of the band. They also have a well defined phase derivative that can be interpreted as instantaneous frequency. This can be useful in a number of applications that require the extraction of such parameters.


Random process $X(t)$ with autocorrelation function given find the mean and the variance Autocorrelation function is
$$R_{xx}(\tau)=\frac{20}{1+2\tau^2}$$
So at $\tau=0$$$R_{xx}(0)=20=E[X(t)X(t)]=E[X^2(t)]$$ The variance is
$$\mathrm{Var}[X(t)]=E[X^2(t)]-E^2[X(t)]=20-E^2[X(t)]$$ As $X(t)$ is WSS (wide sense stationary) the mean is a constant. Is there any way to find its numerical value?

The limit $\lim_{\tau\to\infty} R_x(\tau)$, if it exists, equals $E^2[X(t)]$ and so $E[X(t)]=0$ in this case.

More generally, the mean of a WSS process is nonzero only if the power spectral density has an impulse at the origin. This can be applied to
periodic autocorrelation functions such as $\cos(\omega_0t)$ pointed out
in @MattL's comment. If the Fourier series for a periodic autocorrelation
function has a nonzero DC term, the mean is nonzero.


Generating very short signals So I need to program (using Java) some things that include generated sounds.
I am struggling to understand how to make it work because I require a rather short signal while still using a sampling rate supported by the average sound card.
For example a 100 uSec beep with a 200 uSec delay followed by another beep (initially supposed to be a click but a sine wave also works, would be grateful if you could explain why a sine wave makes it simpler).

My understanding is that I need a high sampling rate (number of samples per second), a frequency within hearing range (for example 1000) and a buffer to hold that data that will be fed to the audio stream provided by the Java libraries.
The issue is that with my limit of sampling rate (hardware limitations) I end up with an almost completely empty buffer because of the tiny intervals of sound it must hold to meet the time requirements.

Am I missing something? Is there another way?
Why does using a sine wave instead of a click make things easier?
Any help would be appreciated.

If you want to generate a very short pulse, it is not usually a problem to have it at any time you wish to have it. That is because modern audio interfaces (sound cards) have an interpolating lowpass filter in the output path, usually with an impulse response that without latency looks like the sinc function:

Figure 1. Sinc function, shown truncated at the ends and critically sampled with the sample number shown in the horizontal axis.

If you play back audio data that is sampled from a time-shifted sinc function with the peak somewhere between sample points, the interpolating lowpass filter will reconstruct the function from the samples so that the peak is in the correct position. Some audio interfaces use a low-latency variant of the lowpass filter that does not have a symmetrical sinc impulse response. Disregarding latency, it has just as good inter-sample timing capabilities if you give it samples of a time-shifted sinc function.

If you need a more rectangular pulse, you could add an analog comparator to the output of the audio interface. You can change the pulse length by varying the voltage level you compare to.

First test if you get that sinc shape from your hardware, like MBaz suggested.


representing simple definite integral as convolution in matlab using CONV A definite integral is a convolution with idial all '1' filter(by the theory)
so i tried to test this issue with MATLAB.
i have taken a function called fun made an integral over it in 1e-6 - 20*e-6 interval the result is called  simbolic_int_result which 37.7
then i have taken fun sampled it in 1e-6 - 20*e-6 interval as shown bellow with 1000 samples and created all "1" filter (idial)
and used CONV to make convolution,but its not working, instead of a scalar i get 2000X1 vector
Where did i go wrong?
Thanks.

c0=2.997*10.^8;% m/s speed of light in vaccum
h=6.625*10.^-34;% J.s Planck constant
k=1.38*10.^-23;%  T/K Boltzmann constant
n=1;

T=307;%  Temperature initiatlisation
fun=@(lambda)(2*pi.*h.*(c0.^2))./((n.^2).*(lambda.^5).*(exp((h.*c0)./(k.*T.*lambda))-1));
simbolic_int_result= integral(fun,1e-6,20e-6)/10;

k=linspace(10.^-6,20*10^-6,1000);
fun_vector=fun(k);

filter=linspace(1,1,1000); % all ones idial filter

convolution_result=conv(filter,fun_vector);

A definite integral is a convolution with idial all '1' filter(by the theory)

No. Where does that theory come from ?

If you want to approximate a definite integral with a discrete sequence, you can simply sum up the samples in the integration interval and scale the sum by the sampling interval.

The output of a convolution of two functions is another function, not a single value.


What is the meaning of negative second for a Morlet wavelet? I kind of familiar with wavelet transform (I just used it in code but do not have a deep understanding about it).
I've been reading an article and came across this sentence:

  where 4Fs is the length of the
  Morlet wavelet with -2 to 2 second support with sampling frequency Fs

and I can't figure out what the -2 second means. By the way the article link is: Deep Convolutional Neural Networks for Brain Computer Interface using Motor Imagery

Continuous wavelets with symmetric envelope are often described, by convention, on a symmetric time interval: $[-T,T]$. The Gaussian being of infinite support, this means it is truncated.

This is only a matter of convention, as the wavelet shape is translated and dilated over the whole time axis. One could thus specify as well the same wavelet on the interval $[0,2T]$.


What is the difference between feature detectors and feature descriptors? What is the difference between feature detectors and feature descriptors? Which among these are detectors and which are descriptors: Harris, SURF, Min Eigen, FAST, SIFT, BRISK

An interest point (key point, salient point) detector is an algorithm that chooses points from an image based on some criterion. Typically, an interest point is a local maximum of some function, such as a "cornerness" metric.

A descriptor is a vector of values, which somehow describes the image patch around an interest point. It could be as simple as the raw pixel values, or it could be more complicated, such as a histogram of gradient orientations.

Together an interest point and its descriptor is usually called a local feature. Local features are used for many computer vision tasks, such as image registration, 3D reconstruction, object detection, and object recognition.

Harris, Min Eigen, and FAST are interest point detectors, or more specifically, corner detectors.

SIFT includes both a detector and a descriptor. The detector is based on the difference-of-Gaussians (DoG), which is an approximation of the Laplacian. The DoG detector detects centers of blob-like structures. The SIFT descriptor is a based on a histogram of gradient orientations.

SURF is meant to be a fast approximation of SIFT.

BRISK, like SIFT and SURF, includes a detector and a descriptor. The detector is a corner detector.  The descriptor is a binary string representing the signs of the difference between certain pairs of pixels around the interest point.


Why do we need motion vectors in MPEG? I'm trying to understand how the encoding of a I and P Frame works.
The I frame is just encoded like a JPEG image.
I'v read that for the P-Frame we need to calculate some motion vectors, apply it to the reference image and calculate the difference between the target and the reference image.
Why do we need such motion vectors? Why can't we just take the difference of the I Frame and the B Frame and encrypt it?
Thanks

Because that's way more efficient at encoding a lot of typical scenes.
Imagine this: you're filming a complex scene with lots of different objects. The camera shifts ten pixels to the right; with your approach, you'll have to save a relatively large difference for every pixel on average, because the pixels 10 px to the left aren't very similar.
When you first encode the fact that this is just a shift by 10 px, that needs nearly no space at all - and then you'll just need to encode how the scene actually changed.