Performing classification based on FFT results

I think the FFT is a bad choice of representation for your problem - it captures many properties of the signal irrelevant to your application, and as you are suspecting, it generates a huge amount of data to process if you extract the FFT of the whole signal.

It seems to me that the most important quantity to consider when studying birdsongs is pitch (fundamental frequency) - all other dimensions of sound (loudness, timbre) are actually variability factors you want to get rid off. For example, two recordings of the same birdsong made in a different environment and with different equipment will exhibit a different frequency response due to the variation in conditions; but fortunately the pitch profile would be exactly the same!

So I suggest you to use a pitch transcription utility (canned solutions : aubio, praat, sonic visualizer...) to extract a pitch contour - a function giving the predominant "note" as a function of time. From that, you could define a feature vector containing pitch statistics (mean, standard deviation, maybe higher order moments); or maybe just build a histogram of pitch values; and this would result in a very compact feature vector suitable for automatic classification. To improve your results, you might then add features capturing the dynamics of pitch over time - dominant modulation rates extracted form the pitch contour, variation of pitch statistics over chunks of a few seconds of audio, etc..


Calculation of an autocorrelation function

I guess you've done almost everything right, but have a problem at the calculation of the expectation value regarding $t$. You should calculate the expectation value of the cosine function. Sadly, it does not simply "go away" as you wrote.

Have a look at the Wikipedia page. There you can find another, more explicit, formula for the auto-correlation function of a function $f(t)$:

$R_{\textrm{ff}}(\tau) = \int\limits_{-\infty}^\infty f(t+\tau)\bar{f}(t) \, \textrm{d}t$ .

(Note, that compared to the Wikipedia page, I have taken the liberty to use the variable $t$ in the integration instead of $u$, which would be the mathematically more accurate version.)

As you can see from this equation, you "integrate away" the dependency on t, and indeed you should be left with a function that is independent of $t$.

Note that there is also an version that does not go to infinite times, but is constrained to a period $T$. Maybe this version is more appropriate in your case. However, the same holds for this version: $t$ is integrated away and should not be a variable in the resulting formula.


How Can JPEG Be Used to Realize a Transfer Application with a Fixed Bit Rate?

You calculated the needed compression ratio correctly.
Now to answer the rest of the question:

The JPEG Algorithm doesn't let you set the size of the output file only the Quality measure. Namely to ensure the size of the output is close to the limit (Namely highest quality) requires trial and error (Though using the Quality you can get a measure of the bits per color).
The JPEG 2000 on the other hand allows you to cut the data at the exact number of bytes required hence makes the task much simpler (You can practically set a threshold on the output size).
The image quality of JPEG at compression ratio of ~100 is really bad (See Wikipedia article on JPEG). While JPEG 2000 can get better results.

To understand the Quality in JPEG, have a look at How JPG Works by Colt McAnlis:

Pay attention to the line:

  Probably most important, is that the quality parameter varies depending on the image. Since each image is unique, and presents different types of visual artifacts, the Q value will be unique as well.

This is what JPEG 2000 makes easier.


Difference Between MFCC of speech and speaker recognition

Both speech recognition and speaker recognition require some set of features to distinguish one speaker (or speech section) from another.

Suppose you have two speakers $S_1$ and $S_2$  and two words $w_1$ and $w_2$. Let's call $u_{ij}$ the utterance of speaker $i$ of word $j$.  Then the MFCC's of the speech are $M(u_{ij})$.

One way to pose the speaker recognition problem is to determine which speaker said a new utterance, $u_{\rm new}$. Assume that the utterance is word $J$. To do this, we just select the speaker such that
$$
I = \arg\min_{\forall i} \lVert u_{iJ} - u_{\rm new}\rVert
$$

How this works will be determined by the inter-speaker variation (i.e. how much the MFCCs of different speakers saying the same word differ).

One way to pose the speech recognition problem is to select:
$$
J = \arg\min_{\forall j} \lVert \bar{u}_{j} - u_{\rm new} \rVert
$$
where $\bar{u}_{j}$ is the average MFCC of all speakers saying word $j$.

How this works will be determined by the inter-word variation (i.e. how much the MFCCs of different words being said by different speakers differ).

So, yes, the MFCCs for both applications are (or can be) the same.


Determining the covariance of point clouds in real-time

To make the Online variance algorithm
work for vectors, change one line:

M2 += np.outer( delta, delta2 )  # deltas 3-vecs, M2 3 x 3

# outer( x, y ) is pairs x[i] * y[j] --
# [[ x0 y0, x0 y1, x0 y2 ],
#  [ x1 y0, x1 y1, x1 y2 ],
#  [ x2 y0, x2 y1, x2 y2 ]]

A small python class for this is under gist.github.com/denis-bz .
It can covariance points with uniform weights (FIR),
or up-weight newer data (IIR), or both.

Is there a tutorial on your application with 3d point clouds ?


Filter banks are one way implement specific discretized versions of a given continuous wavelet transform (or a given SFTF). Their performance depends a lot on how one designs wavelet voices and their subsampling.

if you to the dyadic version (or DWT), then the highest subband spans a $[1/2,1]$ (in reduced frequencies), the next one $[1/4,1/2]$ and so on. So the  frequency resolution is usually poor on the highest subband. Alternatives consist in further dividing each subband (wavelet packets), or in already splitting the data in $M$ subbands instead of $2$ ($M$-band wavelets), or a combination of the two.

Last, the filter banks are quite generic, so nothing prevents your to use an augmented one, with both wavelet-filters and STFT-filers.


80-tap Decimation Filter: divide into stages? / Cost of having a number of taps

Does the number of taps affect a cost of implementation of FIR?

Well, insert your definition of "cost", and your question should really answer itself.

FIR filters belong to the class of linear filters, the combination of N lower order filters can create the desired FIR filter of the higher order. So I was thinking If it makes sense to implement the combination of N lower order filters instead of a one 80 taps filter?

Is the decimation implemented in one stage usually?

Indeed, splitting decimation into multiple stages is something that is commonly done, because then only the first filter needs to run at the highest rate.
Now, 80 taps is not what we'd usually call a large filter, so yes, you can do that splitting, it's even advisable.
However, it might really be premature optimization: Your FPGA is probably fast enough to do all 80 MACs at full rate, and it doesn't even have to do that, seeing that you're building a decimating filter, so everything but the most naive implementation would be a polyphase filter implementation, which reduces your complexity to 80/(decimation) per input sample. Seeing that your filter is only 80 taps long, you're probably also not aiming for a very high decimation, anyway.


Number of microphone limitation for Spherical Harmonics based DOA estimation

First, I do not find any information on what is the order of the sound field. Can anyone elaborate please?

$Y_{lm}(\Omega)$ is the spherical harmonics with mode $m$, order $l$ and frequency $\Omega=(\theta,\phi)$ (dictating angles of arrival, ex: azimuth, spherical, elevation, etc.). According to [1], we have the following formula:

$$Y_{lm}(\Omega) = \sqrt{\frac{2l + 1}{ 4\pi} \frac{(l-m)!}{(l+m)!} }P_l^m(\cos \theta)e^{j m \phi}$$
where $P_l^m(\cos \theta)$ is the associated Legendre functions of degree (mode) $m$ and order $l$. The above serves as a basis of spherical harmonics (similar to Fourier expansion).

  Second, does that mean that I can do very little with 4 microphones? Is there a lower limit on the number of microphones for a DOA estimation and for a number of active sources estimation?

The classical rule of thumb is that given $q$ sources, you need at least $q+1$ antennas to resolve uniquely your sources. However, this is really algorithmic dependent. For example, the paper in [2] (which formulate a 2D AoA DoA estimation) describes shift invariant methods with different lower bounds. So, in principle, the more antennas you have the better resolvability and lower MSE you will have, since the Cramer-Rao bound of DoA estimation is proportional to $\frac{1}{N}$ where $N$ is the number of antennas.

  What other, none neural network, methods can I use to estimate the number of active sources and DOA with (no necessarily by the same model)?

You can use subspace methods like MUSIC, or co-variance fitting. Some generic methods aim at directly optimizing the Maximum likelihood problem like Expectation-Maximization or IQML for specific antenna configurations. A 2D IQML approach was derived in [3].

[1] Rafaely, Boaz. Fundamentals of spherical array processing. Vol. 8. Berlin: Springer, 2015.

[2] Bazzi, Ahmad, Dirk TM Slock, and Lisa Meilhac. "Single snapshot joint estimation of angles and times of arrival: A 2D Matrix Pencil approach." 2016 IEEE International Conference on Communications (ICC). IEEE, 2016.

[3] Bazzi, Ahmad, Dirk TM Slock, and Lisa Meilhac. "Efficient maximum likelihood joint estimation of angles and times of arrival of multiple paths." 2015 IEEE Globecom Workshops (GC Wkshps). IEEE, 2015.


How is a human whistle a pure sine wave?

This is only an approximation!

While the sound of whistling has a very prominent first harmonic, the other harmonics are still present too; there's a non-harmonic noise component; and the first harmonic is more a narrow bump than a sharp line. A better description of the sound of whistling would be white noise filtered through a band-pass filter with a very high Q, and sent through a slight waveshaper adding harmonics. With a slight modulation of amplitude and filter frequency.
A pure sine wave is a  stationary signal, it has constant frequency and constant amplitude. If you observe a whistling sound on a window larger than a few tenths of seconds, you will observe modulations in frequency (because of the melody and articulations between notes), and modulations in amplitude (as you breath and start/end notes).


Maximum cross-correlation coefficient value for time delay estimation

As your plot shows, the second form allows for the correlation peak to be negative. Now, what does a strong negative cross correlation mean? It means the signals are very similar, except one has a negative sign in front of it, i.e., $x_1 \approx -x_2$. Whether or not this makes sense depends a lot on the actual application.

In the application you describe, it would totally make sense to me. See, signals that are transmitted from a source and travel to a destination through some medium are subject to changes in their phase, either due to reflections or the transmission itself. A minus sign corresponds to the opposite phase, which can very likely represent the signal you are looking for, only having gone through an appropriate phase change.

So if I were you, I would allow for the correlation to be negative and look for peaks in the magnitude of the correlation function. Eventually, it still depends on the exact requirements of your application.

If you have the chance, I'd suggest to try some calibration measurements with known source and receiver positions to see what's going on. Then you could even compare both approaches.


Is this jumbled image due to difference between Matlab and Python slicing?

The answer is that the order in which pixels are accessed in Matlab is different from that in Python

The problem is resolved by specifying the order in the numpy reshape function to be the same as that of Matlab using the order='F' flag

Where ‘F’ means to read / write the elements using Fortran-like index order, with the first index changing fastest, and the last index changing slowest

def display_image(input_image) :
  image = Image.fromarray(np.uint8(input_image), 'RGB')
  image.show()

m_resized = np.reshape(np.fromfile('../resnet50_matlab/m_resized.bin', dtype=np.float32),(224,224,3), order='F')

display_image(m_resized)

This results in the correct image being displayed as shown in the image below

Rather than the mangled image displayed using

m_mangled = np.reshape(np.fromfile('../resnet50_matlab/m_resized.bin', dtype=np.float32),(224,224,3))

display_image(m_mangled)


Can we use discrete cosine transform with imaginary part

If we have an OFDM system with $N$ subcarriers, the modulated data based on Fourier exponential function can be expressed as follows:

$S_{ofdm} = \frac{1}{\sqrt T} e^{j2nπFt}$ where $T$ is the OFDM symbol period, $F = 1/T$ denotes the minimum spacing required to guarantee the orthogonality between subcarriers.

Regarding the OFDM-based DCT, the modulated symbol can be done as below:

$S_{DCT-OFDM} = \sqrt {2/T} \cos(2πnft)$

you can compare the two equations. On one hand,  the orthogonality in DCT can be satisfied for a minimum intercarrier distance of $f = 1/(2T)$. On the second hand, the number of subcarriers resulted of DCT are double of those resulted from DFT. So, does that mean we gain or loose half spectrum efficiency?  No, that's not right, because, in DCT, you are just using the real part (or imaginary part) compared with DFT which is using complex basis function. That's all simply how we use DCT with OFDM.

Coming back to your question, is it possible to use the real part and imaginary part?  Yes,  you can use one of them and generally, OFDM based on DCT/DFT should give (almost) similar performance.


Relationship between FFT amplitude and sample size

If you want to have the same amplitudes, you need to scale the FFT result by the length of the FFT:
Fs = 100.0
f = 5;

t1 = np.arange(0, 5, 1/Fs)
t2 = np.arange(0, 50, 1/Fs)

s1 = np.sin(2*np.pi*f*t1)
s2 = np.sin(2*np.pi*f*t2)

plt.figure(figsize=(10,10))
plt.subplot(2,2,1)
plt.plot(t1, s1)

plt.subplot(2,2,2)
plt.plot(t2, s2)

S1 = np.fft.fft(s1) / len(t1)
f1 = np.linspace(0, Fs, len(t1))
S2 = np.fft.fft(s2) / len(t2)
f2 = np.linspace(0, Fs, len(t2))

plt.subplot(2,2,3)
plt.plot(f1, abs(S1))

plt.subplot(2,2,4)
plt.plot(f2, abs(S2))

If you want time and frequency domain signal to have the same power, you need to scale the FFT output by 1/sqrt(length).


What method to use for interpolation and extrapolation?

@Maximilian Matthé, The discrete samples above were made using the function

$$ \mathrm{Truth}(t) = 1.4 + \frac{\sin(13-t)}{13-t} - \sin\big((1.3)(13-t)\big) $$
In this case, the sample rate is several times the Nyquist frequncy of the signal. When that is the case, we can get a very good interpolation over the samples by finding a linear combination of eight sin(f t) terms and eight cos(f t) terms that interpolate the sixteen samples. Results using eight frequencies equally spaced from ($\tfrac{\pi}{16}$ = 0.19635  radians/second) to ($\tfrac{\pi}{2}$ = 1.5708 radians/second) are shown below.

This method is far better for extrapolation (i.e. beyond the samples) than any other method I know of, and the interpolation error (i.e. between samples) is less than $4 \times 10^{-12}$.

This example was not selected to make sure the results look very good, but it is possible to find examples where the approach doesn't excel. If the same frequencies are used and the function is nearly linear over the same sample times, then cubic spline interpolation is better than an interpolating trigonometric series. However, that issue can be overcome by selecting the range of frequencies based on how much oscillation there is between the samples. If the samples are nearly linear, then a set of frequencies closer to zero could be used. Then an interpolating trigonometric series would still work very well if the range of frequencies is appropriate. Also, the samples used in this example are equally spaced, but that is not necessary for an interpolating trigonometric series to work well.


Designing a lowpass filter to minimize aliasing in pre-decimated streaming audio

The best choice of filter depends on your specific application requirements. There are two basic choices: FIR and IIR. IIR will be much more efficient however, it will results in phase distortions. The phase distortions are completely inaudible (unless it's a bizarre outlier case) but clearly measurable. So it depends whether you can tolerate this our not.

In either case you need to decide how close you need to get to the new Nyquist frequency and how much aliasing noise you can tolerate. A typical example would be that you want the passband to extend to 90% of the new Nyquist frequency and that you would like your aliasing products to be below -80dB. Based on these specifications you can then design the appropriate filter. Other considerations include how much pass band ripple you can accept and if you have any constraints on maximum group delay and/or latency.

Here is an example: Let's say you want to downsample from 44.1 kHz to 32 kHz and the new Nyquist frequency is 16kHz. Going to 90% Nyquist (14400 Hz), with 0.1dB pass band ripple and 80 dB of attenuation at 16 kHz could be done with an elliptical filter of 9th order.

As nibot has pointed out, zeroing FFT bins is a poor choice for a low pass filter since the resulting low pass has very big side lobes and aliasing rejection will be quite poor. It would also require a proper implementation of an overlap-add or overlap-save algorithm to deal with a continuous signal.


What sample rate should I use on my ADC?

Depends on your specific application. First of all the guitar has a lot bigger frequency range than you indicate. While the fundamentals are in this range, the guitar has a lot of harmonics and in fact the bulk of the energy is in the harmonics, not the fundamental. So when you plug the low E string at 80 Hz, there is actually not that much energy at 80 Hz but way more at 160Hz, 240Hz, 320 Hz, etc. Same for the higher strings.

In the case it may be simplest to take an off the shelf audio ADC that operates at some standard rate (8000, 11025, 12000, 16000, 22050, 24000, 32000, 44100, 48000). The exact choice  depends on the requirements for your application. For ,say, a tuner, I'd do at least 16k.


The fft provides the spectrum of the audio signal, which tells you the frequencies your signal is made of. In the case of harmonic signals like guitar, violin and other musical instruments, the spectrum exhibits peaks that are (almost) equally separated. For example, 110, 220, 330, 440... Hz. These peaks are called harmonics and they have different amplitudes, which contribute to the perception of timbre.

The perceived note (aka the pitch) depends on the frequency of the first peak, that is the fundamental frequency. It could happen that the fundamental frequency, which is what you want to retrieve with your tuner, is not the strongest one in terms of amplitude. For this reason, you need more complex techniques than just retrieving the maximum of the spectrum. I suggest you to google pitch detection, which is the problem you are working on. There is a lot of material and possible solutions that you can compare and implement.


Octave: Poor performance of convolution by separable 1D kernels

Because loops in MATLAB/Octave are slow (mostly because they are interpreted languages not compiled) and such operations are typically implemented in C/C++ for performance reasons. You could even further speed these up by reverting to better C/C++ implementations and writing your mex wrappers.


Python, resampling tools dont work as expected

The problem is that if your wav file stores samples as integers (typically 16 bit signed), $\texttt{scipy.io.wavfile.read()}$ returns an array of integers. If this is the case, you want to convert your signal into a floating point representation before applying resampling. For a wav file which stores values as N bit signed integers, this can be done by:

// Convert y to float array, with values in the interval [-1.0,1.0]
y_float = float(y / (2.0**(N-1) + 1))


To give a complete answer to this question you're going to need to provide more details about the kind of models you're considering in the first place. But yes, in many cases you can augment those models with a priori constraints on $Y[t]$, such as $0 \leq Y[t] \leq 1$.

For example, if the standard model has some sort of least-squares structure, then adding constraints of that type turns the problem into a bound-constrained least squares problem. There are a variety of approaches to solving such problems, and while they are more expensive than standard least squares, they are quite tractable. And it's very likely that such constraints will produce a better reconstruction.

Even without knowing more, though, I will say this: if your modeling approach does not produce real signals naturally, then you are almost certainly using the wrong modeling approach. It concerns me that you are even proposing taking the real part of the output of some other model. You should be searching the space of real signals if you know that to be the underlying structure.


harmonic product spectrum on guitar string signals

HPS (harmonic product spectrum), ASDF, AMDF, and autocorrelation are all similar periodicity estimation methods (they use different weightings and algorithms), which might be useful for guitar note pitch estimation.  These methods have different estimation accuracy/error statistics, which may vary depending on the types of guitars, strings, microphone setup, and etc., tested.  Experiment.

HPS has one interesting advantage.  The overtone series of the lower guitar strings can be slightly inharmonic, partially due to the physics of actual guitar strings, which have non-zero thickness and stiffness.  However, there is some experimental evidence that the stretch between the inharmonic overtone series might actually correspond more closely to the typical human psychoacoustic perception of pitch, than the actual fundamental vibration mode frequency of the guitar string.

That said, naive use of these methods may be too heavyweight for an Arduino AVR processor.  Perhaps an Arduino using an ARM M0 would be more suitable.  Also, to reduce the computational load of doing a full FFT or full autocorrelation-equivalent, you may want to use a modified zero-crossing pitch estimator.  Just make sure to use a large enough set of zero crossing span pairs (not just adjacent crossings) so that you can use something like a single AMDF to pick which set of span pairs represent overtone, harmonic or noise lags, and which span pairs might represent fundamental pitch lags.  Once you validate your set of span pairs or lags, you may want to average, or median filter.


Matlab: Different 1-D FFT results with and without zero-padded

I don't know how MATLAB handles this internally. What will help in your situation however is multiplying your unpadded signal by a suitable window function, like a hann window.

What you see in your diagrams is a convolution with a $\text{si}$ function, which is the fourier transform of a rectangle function. If you just take the signal as it is and transform it, you essentially multiply it by a rectangle function.

If you try this, for both of your cases, you should get matching results.


azimuth extent for range migration algorithm

I'm going to assume that you are operating in Stripmap mode with no squinting i.e. broadside mode. The aperture size is given by the real antenna azimuth beamwidth (in radians) times the range to the target. What you use as the beamwidth is really a design question - you can use 3dB point or the null-null width. The beamwidth will also be affected by the elevation beampattern somewhat.

You can see near range targets will have a smaller synthetic aperture and thus have smaller coherent gain ( fewer azimuth pulses that can be coherently added). Farther range targets will have a larger coherent gain because more pulses can be coherently added. Most SAR systems have a relatively small range swath widths and so the change in the aperture length is pretty small. The beamwidth is what limits the extent of the azimuth Doppler freqency and thus determines the azimuth resolution.

There are a couple of approaches to handling larger range swath widths. The problem with large range swath widths has more to do with handling the change in range curvatures over the range swath rather than the change in the aperture length. One of the better books to discuss this is by Cumming and Wong:"Digital Processing of Synthetic Aperture Radar Data"

Note that in practical processing systems - if the synthetic aperture at your range of interest contains N pulse, you tend to process much more than N pulses so that you get more azimuth data at the output of your processing. If you only process N pulses - then you get 1 point in azimuth, while if you process 2N pulses you can process N+1 synthetic apertures (1->N,2->N+1, ...,N+1->2N) and you get more effective processing.

In Spotlight mode you steer the antenna beam at one spot on the ground throughout the data collection. This allows a higher-resolution image to be formed but over a smaller area.


Couple of approaches come to mind: you could pick an arbitrary "first" red pixel, then take a square of say 5×5 pixels around it, and simply figure out which quadratic function is a best fit for these, and infer the curvature from that.
Pick an arbitrary direction to start from that point and pick the next pixels; calculate the curvature of the 5×5 square surrounding that pixel. Find the absolute difference between last and current curvature; if it changed too much, you've found an abrupt corner.
In your specific example, applying to the whole image a 5×5 moving average filter would also suffice: Far as I can tell, there shouldn't be more than a relatively low number of red pixels in every 5×5 square, so if there's one with more, you've got your type of defect.
Other options include morphological classical tricks like erosion / dilation: use one operation to surround every red pixels with a fixed-diameter disk of red. Then do the same, with a diameter one linewidth (or more) larger disk, with the white in the picture. Check whether anything remains red afterwards.
PS: JPEG is the wrong image format / compression for line-style graphics. You'd want to work with PNG.


Visualisation of bandwith of filter with sampling

$B$ is the cut-off frequency of the anti-alias filter needed before sampling; $\omega_0$ is the signal bandwidth, and $\omega_s$ is the sampling frequency.

$B$ can also be the cut-off frequency of the (ideal) reconstruction filter (band-limited interpolator).

If this is the case, this suggests that the reconstruction filter is an ideal low-pass filter with cut-off frequency between the signal bandwidth and the next overlapped frequency component ($\omega_s - \omega_0$). The sampled signal contains repetitions of the original specturm every $\omega_s$, so when reconstructing, you want to get rid of all these repetitions and keep only the spectrum centered at 0.

I don't know what literature you are referring to, and why is it so crappy.


Why does each node in the wavelet scattering transform split into multiple paths?

We take a wavelet transform over other wavelet transforms, which involves multiple wavelets. E.g.:

First-order uses 100 wavelets to tile the frequency domain and produce 100 coefficients
Second-order uses 16 wavelets to tile the frequency domain and produces 8 coefficients per first-order coefficient, yielding a total of 800
Third-order uses 9 wavelets ... etc

(Note second-order wouldn't actually produce 800 in implementation since xi2 > xi1 wavelets are discarded, where xi == center frequency).
The idea is, first-order coefficients capture temporal modulations (AM-FM). Second-order capture the modulations of those modulations, extracting higher order structure, and the process repeats. Progressively we require less wavelets as there's less fine-grained structure to capture (also energy and non-redundancy reasons), but it's still more than one wavelet (i.e. path), as just one wavelet cannot cover entire frequency domain without losing information (unless it's very wide, but that's a bad wavelet).
See helpful lecture, and Python implementation. I also might write a detailed post on this topic in coming weeks.


Is there a way to break an analog sound signal up and send it through two radio chips at double the data rate of one?

This is problematic because I don't see any guarantees that a pair of those chips will provide you with same end-to-end latencies. You can probably cheat a bit if you clock both chips in the pair with the same clock, and provide the pair with a common synchronous reset signal.

Assuming the end-to-end latency is equal down to a single sample, you can leverage the fact that this chip family can take digital audio input and provide digital audio output. Thus you can sample the audio at a double data rate, and split it digitally between two chips using a bit of discrete logic or a small CPLD. The reverse of this process can be applied at the receiving end.

If you wish to align the signals at the receiving end in spite of no provisions for this by the chipset itself, you can be inserting synchronization marks (single-sample clicks) at the transmitting end. Those clicks are easy to detect on the receiving end, and can be "fixed" by interpolating over two adjacent samples. If the clicks are sent on a random schedule - say uniformly distributed between 0.5 and 2 seconds, they won't appreciably affect the signal. Any errors from the click removal will be spread-spectrum.


You can't identify aliasing with a simple sinusoid at a specific frequency and that in a way is the whole point about trying to avoid it. You can't know if the sinusoid you are "looking at" is $Q$ Hz or $2Fs-Q$Hz.

A single aliased sinusoidal component looks just like a non-aliased sinusoid. If you want to experience aliasing, you have to attempt it with either a more complex waveform or a sinusoid that is changing in time.

One way to "experience aliasing" is to undersample a chirp in the following way:

Fs = 8000;t=0:(1./Fs):(5-1./Fs);p=2.*pi.*t; %Sampling Frequency, Time Vector, Phase Vector
y1 = chirp(t,0,5,Fs/2.0); %Create a chirp that goes from DC to Fs/2.0
spectrogram(y1); %Have a look at it through spectrogram, please pay attention at the axis labels. This is basically going to be a "line" increasing with time.
soundsc(y1,Fs); %Listen to it...It clearly "goes up" in frequency
y2 = chirp(t,0,5,Fs); %Now create a chirp that goes from DC to Fs
spectrogram(y2); %Have a look at it through spectrogram
soundsc(y2,Fs); %Listen to it...Do you "get" the folding of the spectrum?

In general, you can think of Sampling as Modulation because this is what is effectively happening at the sample and hold element of an ADC converter.

This will enable you to understand more easily concepts like undersampling for example (and applications where it is perfectly OK to sample at lower than the Nyquist rates). But also, you can load a WAV file in MATLAB (with 'wavread') containing some more complex signal and prior to listening to it with 'soundsc', simply multiply it with a "square" wave* (you might want to lookup function 'square') at some frequency lower than the WAV file's $F_s$. This will effectively introduce the key (unwanted) characteristic of aliasing which is this folding of the spectrum. The result is not very pleasant so you might want to keep your speakers volume down.

I hope this helps.

*EDIT: Obviously, "square" returns a square with an amplitude in the interval [-1,1] so prior to multiplying it with your signal it would be better to rescale it as:

aSquareWave = (square(100.*p)+1.0)/2.0 % Where p is the phase vector and here we are producing a square wave at 100Hz (given an Fs of 8kHz as above). aSquareWave's amplitude is now in the interval [0,1]


How to detect if a signal raises sharply or slowly?

I believe that you can calculate the derivate which is the rate of signal change. In discrete signals this is obtained (approximated) by:

$$x'[k] = x[k]-x[k-1]$$

(by other words, it's the difference between consecutive samples)


Sampling a signal with varying frequency

If your frequency increase linearly, then you have a chirp signal with linearly increasing frequency.

https://en.wikipedia.org/wiki/Chirp_spectrum

There is a closed-form solution for the Fourier transform. I'm not sure how that will help you. Why not pick a sampling frequency

$\omega_s \gt 2 \omega_{max}$

Since your highest frequency is $\omega_{max}$, no point in considering lower frequencies for your sampling rate. Is this for a control-loop algorithm? In that case, I would use a more rule-of-thumb criterion.

$\omega_s \gt 20 \omega_{max}$

This makes you less susceptible to sampling delays (among other things)


Let me provide a method, applied only for finding the impulse response $h(t)$ of an LTI system characterised by an LCCDE of the form $ \sum_{k=0}^{N}{ a_k  {{d^k y(t)}\over {dt^k}}} = \sum_{k=0}^{M}{ b_k  {{d^k x(t)}\over {dt^k}}}$ by using the classical time domain approach. A method which is ignored as it's replaced by transform domain techniques instead...

Note that the input $x(t) = \delta(t)$ is formally a problematic signal (function) and therefore mathematically oriented classical books on ODEs, do generally avoid any discussion of such functions and their solutions, unless the scope of the book specifically includes generalised functions (distributions), usually to be used in some engineering or physical fields.

The procedure: Consider an LTI system which is causal with initial rest conditions. We say that when the input to this system is $x(t)=\delta(t)$ an impulse, then its output $y(t)=h(t)$ is the impulse-response of the system. We find this solution $h(t)$ in two steps, by breaking the LCCDE into two parts as inspired by a serial (cascade) implementation of two LTI systems corresponding to the right and left sides of the LCCDE.

I/O relationship of Part-I is given by the LCCDE: $$\sum_{k=0}^{N}{ a_k  {{d^k y(t)}\over {dt^k}}} = x(t)$$ which represents the first stage with its input $x(t)=\delta(t)$ and we denote its solution (stage1-output) as $h_0(t)$ which is actually the impulse response of Part-I.

And the I/O relationship of the Part-II is given by the equation: $$ y(t) = \sum_{k=0}^{M}{b_k{{d^k x(t)}\over {dt^k}}} $$ which requires nothing but simple summation of its input, $x(t) = h_0(t)$ ,and its derivatives to compute the output as $$h(t) = \sum_{k=0}^{M}{b_k {{d^k h_0(t)}\over {dt^k}}}$$ Therefore we need to find $h_0(t)$ of the Part-I to simply compute the impulse response $h(t)$ of the complete system.

To find the solution of part-I: $\sum_{k=0}^{N}{ a_k  {{d^k y(t)}\over {dt^k}}} = x(t)$ whose solution $y(t)$ is what we call $h_0(t)$ , we we make the following observation that $y(t)$ is composed of two parts as $$y(t) = y_{h}(t) + y_{p}(t)$$ where $y_{h}(t)$ is the homogeneous solution corresponding to $x(t)=0$ for all t, and $y_{p}(t)$ is the particular solution corresponding to input $x(t)$ for $ t > 0$, in particular as what remains as input when time goes to infinity.

Now, we can see that for the input $x(t)=\delta(t)$ the particular solution $y_p(t)$ for $t>0$ is zero as the input is zero for $t > 0$. Then the output simplifies to $h_0(t) = y(t) = y_h(t)$

Now, we have to find the homogeneous part of $y(t)$ as the solution of the equation $$\sum_{k=0}^{N}{ a_k {{d^k y(t)}\over {dt^k}}} = 0$$
following the usual procedure by finding the roots of its characteristic polynomial in complex-s: $$p(s) = \sum_{k=0}^{N}{a_k s^k}$$

which is obtained when we simply insert $y_h(t) = Ke^{st}$ as a solution canditade to the homogeneous differential equation $\sum_{k=0}^{N}{ a_k {{d^k y(t)}\over {dt^k}}} = 0$. The solution is then hypothesised to be $$y_h(t) = \sum_{k=0}^{N}{ A_k e^{st}}$$ for the case when all the roots are distinct. If they are not, you know what to do right ;)
The solution therefore now requires to determine the N unknown coefficients $A_k$ from N-inital conditions, at time t=0 for the framework of this method, possesed by the output $y(t)$ and its derivates up to order N-1.

Remember that this system was defined to be under inital rest with complete zero inital conditions prior to the application of any excitation. Hence if we were to solve the homogeneous solution $y_h(t)$ for an ordinary (formally valid) input signal $x(t)$ to the system, the homogeneous solution will be identically zero, as the initial conditions would all be. (all the coefficients $A_k$ of the exponential terms turn out to be zero). But the magic of the generalised input $\delta(t)$ is that it will set one of those inital conditions to non-zero, thereby, enabling a non-zero homogeneous response to exist, which will become the solution of the part-I as well.

It can be shown that the inital conditions of the part-I due to the excitation $\delta(t)$ will be set as $y(t)=0 ~,~ y(t)' = 0 ~,~ y(t)''=0 ~,~... ~,~ y(t)^{N-1}=1/a_N$ all at $t=0$.

From this information you can proceed to solve the algebraic equaitons as usual and find $h_0(t)$ from which you can compute the impulse response $h(t)$ for all $t > 0$. What happens for all $t <0 $? As we stated, this system was LTI and causal with initial rest condition. For such a system if the input $x(t)$ is zero for all $t<0$ then that output is also zero for all $t<0$ which is achieved by multiplying $h_0(t)$ via $u(t)$ to represent the solution of Part-I for all t.

It is instructive to apply this procedure to your example:
$$y'' + 5 y' + 6y = -x'$$

Lets first solve Part-I:  $$y'' + 5 y' + 6y = x$$
To find $h_0(t)$ when the input is $x(t)$ = $\delta(t)$ first solve the characteristic equation: $s^2 + 5s + 6 = 0$ which yields $s=-2$ and $s=-3$.

The proposed homogeneous solution is therefore $$y_h(t) = A_1 e^{-2t} + A_2 e^{-3t}$$

For which we need 2 inital conditions to find $A_1$ and $A_2$.
These initial conditions are set by the impulse $\delta(t)$ to be: $y(0)=0$ and $y(0)'=1/1=1$. Using these initial conditions we find $A_1$ and $A_2$ to be $A_1=1$ and $A_2=-1$ and this yields the solution to part-I as
$$ h_0(t) = [ e^{-2t} - e^{-3t}] $$ for all $t > 0$.

Also note that this system was LTI and causal with initial rest and input $x(t)=\delta(t)$ is zero for all $t<0$ implies that output is also zero for all $t<0$ which is simplified by using a unit step $u(t)$ function to format the output as:
$$ h_0(t) = [ e^{-2t} - e^{-3t}]u(t)$$ for all t.

Finally, to solve part-II and to get the impulse response of the overall system, apply this $h_0(t)$ as an input to Part-II whose output is simply

$$y(t)= -x(t)' \longrightarrow h(t)= -h_0(t)'$$

$$h(t) = -[ -2 e^{-2t} + 3 e^{-3t}]u(t) - [e^{-2t} - e^{-3t}]\delta(t) $$
$$h(t) = -[ -2 e^{-2t} + 3 e^{-3t}]u(t) - [1 - 1]\delta(t) $$
$$h(t) = -[ -2 e^{-2t} + 3 e^{-3t}]u(t) - 0\delta(t) $$
$$h(t) = [2 e^{-2t} - 3 e^{-3t}]u(t)$$

Note that the impulse at the output is cancelled by the property $f(t) \delta(t-a) = f(a)\delta(t-a)$


At a high level, it's a relatively simple tradeoff: software development for DSP or general-purpose microprocessors is typically much less labor-intensive than FPGA development. At the same time, the amount of computational throughput offered by an FPGA typically far exceeds what a CPU can do, because you can tailor the logic specifically for your application.

Most any processor that is advertised for "DSP" will have multiply-accumulate (MAC) instructions available, and it is common for DSP processors to be able to execute multiple instructions per cycle. If you read the specifications for TI's C6713, you'll see that they advertise that it can complete up to 8 instructions per cycle (not necessarily all MACs, but a mix of several types of instructions). Note that it is pretty difficult to get close to that throughput, but it is possible if you invest enough effort learning how the architecture works and structuring your code accordingly.

One other major difference between an FPGA implementation and the particular DSP that you noted is that the C6713 supports native floating-point capability (instead of just fixed-point). For some applications, this can make algorithm implementation much simpler, as you don't have to worry about the issues that come along with fixed-point computations. While it is possible to implement floating-point computations in an FPGA, such an implementation would be much more demanding on the hardware, so depending on how complex your algorithm is, you might run into issues with resource utilization, achievable clock speed, and so on.

You should start by quantifying what kind of computational demands your algorithm is going to place on your hardware. Then, decide how much engineering time you can allocate to the project; from there, the choice will likely be clear.


This is the classical use case for a Phase-Locked Loop (PLL). You can control speed vs noise robustness by adjusting the loop filter bandwidth.


Can I compute signal-to-noise ratio on a vector-valued or multi-channel signal?

You can take the way SNR is calculated in RGB images as a starting point.

In RGB images, the image SNR is calculated using the familiar formula

$$
\frac{\mu}{\sigma}
$$

Where $\mu$ is the "whole image" mean and $\sigma$ is the "whole image" standard deviation. To calculate these "whole image" values, each of the RGB channels is first converted to the common currency of luminosity by one of several formulas, generally weighted sums:
$$
L = W_{R}\mu_{R} + W_{G}\mu_{G} + W_{B}\mu_{G}
$$

Where $\mu$ are means and $W$ are respective weights for the R, G, and B channels of the image (in the Adobe formula, for example, weights are 0.2974, 0.6273, and 0.0753 respectively), and L is luminosity.

For whole image standard deviation, the variance includes the covariances of the channels, so the full formula is
$$
Var(L) = \sum_{i,j} (r_{i,j}^2Var_{R} + g_{i,j}^2Var_{G} + ...
$$
$$
b_{i,j}^2Var_{B}+ 2r_{i,j}g_{i,j}Cov_{R, G} + ...
$$
$$
2r_{i,j}b_{i,j}Cov_{R, B} + 2g_{i,j}b_{i,j}Cov_{G, B} )
$$
where $i$ and $j$ are pixel indices, which is a more expensive calculation and not as commonly used in industry. Oscar de Lama reports that the Adobe formula simply disregards the covariances:
$$
W_{R}^2Var_{R} + W_{G}^2Var_{G} + W_{B}^2Var_{B}
$$
and suggests a compromise formula:
$$
W_{R}Var_{R} + W_{G}Var_{G} + W_{B}Var_{B} ...
$$
$$
+ (2/3)(W_{R}\mu_{R}^2 + W_{G}\mu_{G}^2 + W_{B}\mu_{B}^2) ...
$$
$$
- (2/3)⋅(\hat{W}_{RG}⋅\mu{R}⋅\mu{G} + \hat{W}_{RG}⋅\mu{R}⋅\mu{B} + \hat{W}_{RB}⋅\mu{G}⋅\mu{B})
$$
Where $W$ are the earlier weights, and $\hat{W}(X,Y) = \sqrt{XY} / K$ where $K$ normalizes $\hat{W}$ such that $\sum \hat{W} = 1$.

If the vector-valued image weights all components equally, then W reduces to 1, simplifying the formulas.


Is Spectral Coherence (MSC) dependent to the initial amplitude of the comparing signals?

From the definition of the (magnitude-squared) coherence

$$C_{xy}(f)=\frac{|G_{xy}(f)|^2}{G_{xx}(f)G_{yy}(f)}$$

with the cross-spectral density $G_{xy}(f)$, and the power spectra $G_{xx}(f)$ and $G_{yy}(f)$, respectively, it is clear that scaling of $x(t)$ or $y(t)$ does not change the value of $C_{xy}(f)$, because the scaling constants appear in the numerator as well as in the denominator and cancel out.

The coherence is a normalized quantity (and hence scaling-invariant), and it always satisfies $0\le C_{xy}(f)\le 1$.


Can STFT (Short-time Fourier Transform ) be more useful than FFT for analyzing stationary signals under some circumstances?

In the continuous version of the short-time Fourier transform, if you choose a uniform unit-height window, it just boils down to the continuous Fourier transform.

In practice, for the discrete case, you can find STFT parameters that just give you the same information as that of the FFT. Hence, STFT can always be as useful as the FFT.

Since real signals are often finite length, they are not stationary in the purest sense. Even if the noise is "stationary", we only have access to few realizations.

For instance, the estimation of the parameter of an additive noise from a few realizations would probably be more robust in the STFT space than from the FFT.


Why are LTI systems defined by convolution, why not in any other way?

it takes more work to do this for continuous-time systems, so i'll just do it for discrete-time and we can all wave our arms to extend to continuous-time.

a general discrete-time input signal can be expressed as

$$ x[n] = \sum\limits_{i=-\infty}^{\infty} x[i] \delta[n-i] $$

where $\delta[n]$ is the Kronecker delta or "unit impulse function"

$$ \delta[n] \triangleq
\begin{cases}
1, & \text{if }n = 0 \\
0, & \text{if }n \ne 0
\end{cases} $$

so think of your input $x[n]$ as the sum of a bunch of individual impulse functions, each delayed to time $i$ and weighted by a constant which is equal to $x[i]$.

so, let's say we have a Linear and Time-Invariant discrete-time system, represented with operator, $\mathbb{T}$, and output:

$$ y[n] = \mathbb{T}\left\{ x[n] \right\} $$

and we input into that system a single unit impulse, $\delta[n]$.  by definition, the "impulse response" of the LTI system is the output coming out when a unit impulse, $\delta[n]$, is input.

$$ h[n] \triangleq \mathbb{T}\left\{ \delta[n] \right\} $$

now, imagine a general input, $x[n]$ is input to $\mathbb{T}\{ \}$.  then, because it's Linear:

$$ \begin{align}
y[n] & = \mathbb{T}\left\{ x[n] \right\} \\
& = \mathbb{T}\left\{ \sum\limits_{i=-\infty}^{\infty} x[i] \delta[n-i] \right\} \\
& = \sum\limits_{i=-\infty}^{\infty} \mathbb{T}\left\{ x[i] \delta[n-i] \right\} \\
& = \sum\limits_{i=-\infty}^{\infty} x[i]  \mathbb{T}\left\{ \delta[n-i] \right\} \\
\end{align} $$

and because it's Time Invariant,

$$ h[n-i] = \mathbb{T}\left\{ \delta[n-i] \right\} $$

then, putting the two together, you get

$$ \begin{align}
y[n] & = \sum\limits_{i=-\infty}^{\infty} x[i]  \mathbb{T}\left\{ \delta[n-i] \right\} \\
& = \sum\limits_{i=-\infty}^{\infty} x[i]  h[n-i] \\
\end{align} $$


Reason for bimodal behavior while low second fourier coefficient

Two peaks does not imply a frequency of 2 - for example the peaks could be very narrow and/or their spacing could corresponding to a much smaller interval than half the time-series length.

There is no information about the overall shape of the time-series that can be drawn from the knowledge of the second DFT coefficient.

I have plotted here 15 time series (left) and the magnitude of their FFT (right). 5 unimodal, 5 bimodal, and 5 trimodal series.


By definition:
$$x(t)*h(t)=\int_{-\infty}^\infty x(\tau)h(t-\tau) \ \mathrm{d}\tau$$
If any of the functions inside the integral equals 0, then the product does as well. We know that:
$$x(\tau)\neq0  \iff \tau\in(0,1)$$
$$h(t-\tau)\neq0  \iff t-\tau\in(0,a) \iff \tau \in (t-a,t)$$
So the convolution integral can be simplified:
$$x(t)*h(t)=\int_{-\infty}^\infty x(\tau)h(t-\tau) \ \mathrm{d}\tau
= \int_{\max(0, t-a)}^{\min(1,t)} x(\tau)h(t-\tau) \ \mathrm{d}\tau
$$
Just as you had found out.
Notice that this works just because $0<a\leq 1$. If that was not the case, then some other things should be taken into account.


If the input sequence is band limited to $\pi/2$, the second alternative is also OK (as you realized with your argument). However, if it not band limited you will get aliasing, so it is not a feasible solution in the general case.


How to remove text region in a document image? or extract image from a document image

What you are asking is called 'page layout analysis' in OCR. Take a look at https://en.wikipedia.org/wiki/Document_layout_analysis.


Where does the delta function come from if we derive autocorrelation directly?

The book is indeed not consistent. Continuous-time (zero-mean) white Gaussian noise has an autocorrelation function

$$R(\tau)=\sigma^2\delta(\tau)\tag{1}$$

and a constant power spectral density (PSD)

$$S(\omega)=\mathcal{F}\{R(\tau)\}=\sigma^2\tag{2}$$

The power (or variance) of such a process is

$$R(0)=\int_{-\infty}^{\infty}S(\omega)d\omega=\infty\tag{3}$$

So in the book the variable $\sigma^2$ denotes two different things: the finite value of the (constant) PSD, and the (infinite) variance of the process. But it can't be both at the same time. So the equation $R(0)=\sigma^2$ doesn't make sense, if at the same time we have $R(\tau)=\sigma^2\delta(\tau)$.

EDIT: In reaction to your comment concerning the definition of the auto-correlation function: the definition given in the book is correct (for real-valued wide-sense stationary processes):

$$R(\tau)=E[x(t)x(t-\tau)]\tag{4}$$

Since we know that $x(t)$ is white and zero-mean, we know that

$$R(\tau)=E[x(t)x(t-\tau)]=0,\quad\tau\neq 0\tag{5}$$

For $\tau=0$ we get

$$R(0)=E[x^2(t)]\tag{6}$$

which is the power of $x(t)$, and, since $E[x(t)]=0$, it also equals its variance. And now comes the error in the book: the authors use $\sigma^2$ to denote the variance: $$R(0)=\sigma^2\tag{7}$$ But then they claim

  The autocorrelation function of white noise must therefore be a delta function. $$R(\tau)=\sigma^2\delta(\tau)\tag{8}$$

But this is inconsistent with Eq. $(7)$ (of this answer). Note that both $(7)$ and $(8)$ are correct, but together they are wrong because they use the same symbol $\sigma^2$ with a different meaning. In $(7)$, $\sigma^2$ is the variance of $x(t)$, whereas in $(8)$, $\sigma^2$ is the value of the (constant) PDS of $x(t)$.

In the given example, $x(t)$ is white, so its power (variance) is infinite. This means that the constant $\sigma^2$ in $(7)$ cannot be finite. On the other hand, the constant $\sigma^2$ in ($8$) is finite, because it is simply the value of the constant PSD of $x(t)$.

The fact that the auto-correlation function of a zero-mean white wide-sense stationary process must be a weighted delta impulse at $\tau=0$ is obvious from the requirement that the PSD must be constant. The auto-correlation is the inverse Fourier transform of the PSD, and the inverse Fourier transform of a constant is a Dirac delta impulse.


Periodicity of the output of causal time-invariant systems when the input is periodic

Your line of reasoning is correct. For a (not necessarily linear) time-invariant system the following always holds: if $y(t)$ is the response to an input signal $x(t)$, then $y(t+T)$ must be the response to the shifted input signal $x(t+T)$. If the input is periodic with period $T$, this implies that $x(t)=x(t+T)$, and, consequently, $y(t)=y(t+T)$, i.e. the output signal is periodic with the same period as the input signal.

Furthermore, if the system is causal and the input up to time $t_0$ is periodic, the output up to time $t_0$ must be periodic too, because a causal system cannot know the (possibly non-periodic) input beyond time $t_0$.

This is (part of the) theory. But then comes @MBaz with his "counter-example" of an ideal integrator. The problem here is indeed, as already noted in a comment, the integrator's instability. Showing the time-invariance of a linear system involves manipulating the convolution integral. However, if for the given input signal that integral is divergent, the proof becomes invalid. And this is the case with systems having poles on the $j\omega$-axis (or, for discrete-time systems, on the unit circle). If an LTI system has a pole on the $j\omega$-axis and if you excite it with the corresponding pole frequency, you get an output signal with increasing amplitude, which is obviously non-periodic. You can also see it in the following way: the response of an LTI system to a complex exponential $e^{j\omega_0 t}$ is

$$y(t)=H(\omega_0)e^{j\omega_0 t}\tag{1}$$

Eq. (1) of course requires that $H(\omega_0)$ can be evaluated. For an ideal integrator this is not the case for $\omega_0=0$. Note that if the input signal to the integrator had no DC component, then the response to a periodic input signal would indeed be periodic.


Spectrogram from Complex Morlet wavelet ( recreating article result from data)

Well I managed to find this MATLAB package

https://sites.google.com/site/rwfwuwx/Home/mfeeg

It does almost the same operation that as been done in the paper..

the related files in the package are:
mf_tfcm.m & mf_cmorlet.m


Pole-Zero Configuration of Glottal Signal

Recall what a pole is? It's where the z-transform tends to infinity. What value(s) of $z$ would cause $G(z) =\frac{1}{(1-\beta z)^2}\rightarrow \infty$? (This will answer your reciprocal question).

Next, note that $\beta <1.$ So where would $\beta$ go relative to 1?
(This will answer your question about why the pole is drawn in the location shown.)


Sound glitches using fft segmented convolution

One possible way to reduce cracking noises between convolutional block filtering parameter changes is to use an overlapping of the blocks and do a cross fading (linear or raised-cosine) of the overlapping segments.


Is it always necessary to use two cameras to find the distance between camera and object?

There are various ways to measure distance and 3D using a single camera.
These include:

Moving the camera and taking another image - this simulates a 2 camera system.
If you are imaging an object of a known size, then you can know the distance with some simple trigonometry.
You can project a known pattern from a known projector location to reconstruct 3D. This is structured light and e.g. how Kinect works. Obviously the image(s) with the pattern will appear different from the normal image, so multiple exposures might be neccesary if you also want the appearance image.
You can use Depth-from-Defocus by capturing multiple images at multiple focii.
Sometimes an auto-focus lens can provide rough estimates on the depth of the object in focus. If the object is big enough in the frame this may be sufficient for some applications.

Each method has different trade-offs. You must consider which one suits your app.


Python: How to calculate image perimeter orientation histogram?

I would like to find the perimeter orientation histogram of an object inside an image.
  The object is returned into an array but I don`t know how to find the orientation of its perimeter.

As far as the perimeter edge orientation histogram is concerned, this is a histogram of a set of integer numbers in the range $[0 \ldots 5]$. In section 3.3 it is explained that the orientation of the perimeter "edges" is quantised in 6 directions, each one represented by oriented Gaussian Filters. The orientation of the edge is taken as the "strongest" output of one of these filters.

  That is to return the perimeter pixels coordinates as a 2-D array of X and Y. What I found in Python is to just return the total number of pixels in the object perimeter but no information about the coordinates of perimeter.

This has nothing to do with Python. Finding the perimeter is not a matter of "calling a function". It is a set of standard techniques and the key towards them is the phrase "Perimeter - total number of pixels on the segment perimeter, computed with morphological operations", again, from paragraph 3.3 of the paper cited.

  Can anyone help me find the perimeter coordinates of an object and then calculate the orientation of the perimeter pixels?

You first need to understand morphological operators. This is a great resource and so is this one, particularly section 3.6 which is exactly what you are trying to do here.

In brief, given a thresholded binary image where $1$ denotes a pixel in the Foreground and $0$ denotes a pixel in the Background, the perimeter of it is defined as the set of pixels that connect the Foreground with the Background regions. Obviously, here, Foreground, Background, 1,0 are interchangeable and depending on context. But for the purposes of this discussion we define Foreground as the pixels that the segmentation algorithm has determined to belong to some identifiable area, Background as the pixels that the segmentation algorithm has determined DO NOT belong to some identifiable area and 0,1 are just two numbers we use in place of Foreground and Background (It could be a different integer for each identifiable segmentation region).

The brute-force way to do this is to scan through the "patch", examining the 4 or 8 connectivity of each pixel and looking for pixels that are not fully surrounded by Foreground pixels. THOSE PIXELS are the ones that exist in the perimeter of a region identified by the segmentation algorithm.

OK, so the next question here is Where do i get these regions identified? and the answer is, from the dynamic texture modelling step that is described in section 3.1.

What is dynamic texture modelling?

It is a technique by which an image is modelled as a set of particles that tend to move in well defined patterns. From this perspective, a dynamic texture model is trying to discover the pattern of movement that these particles follow. Having discovered this, it is possible to "judge", if some particle's movement belongs or not to a specific model.

In other words, when you observe scenes of cars moving on a motorway, birds flying in the sky, particles flowing in some liquid, or even...pedestrians moving along predetermined paths through a park, then these particles tend to follow well defined routes. For more information, please see this piece of work

This is how the authors are discovering those grayscale blobs on Figure 3.

Specifically:

"We adopt the mixture of dynamic textures [18] to segment
the crowds moving in different directions. The video
is represented as collection of spatio-temporal patches (7 ×
7 × 20 patches in all experiments reported in the paper),
which are modeled as independent samples from a mixture
of dynamic texture models [19]. The mixture model is
learned with the expectation-maximization (EM) algorithm
[18]. Video locations are then scanned sequentially, a patch
is extracted at each location, and assigned to the mixture
component of largest posterior probability."

So, having:

Segmented the crowd on each frame by first having learned the pattern of movement of the pedestrians.
Discovered the pixels that lie on the perimeter of those patches

You then take the position of each pixel on the perimeter of the region, center six differently oriented Gaussians on top of it and select the "strongest" direction from these six.

In the end of this process, you will have the set of directions required to construct the histogram.

In terms of Python implementations of these techniques, please see this link and this link.

Hope this helps.


redundancy of sin and cos waves with real data

Both sin and cos waves are needed. Sure, you can find a signal (somewhere between a sine and a cosine, which seems as if it can be evaluated using either a sine or a cosine but in the general case you need both waves. And FT needs both waves anyway because internally it works with complex numbers regardless of the input  you pass to it (real input is usually interpreted as real0, img0, real1, img1...)

Suppose you have a pure sine wave and want to estimate its spectral content. Which wave would you use, a sine or a cosine?

It should be evident that only sine waves can accurately estimate the frequency and amplitude of your tested sine signal (provided that they're in phase and the frequency resolution is narrow enough to reduce spectral leakage effetcs). If you now inspect the DFT/FFT results, you'll see that only a single, imaginary coefficient (corresponding to the sine wave being tested) contributes to the overall amplitude/magnitude.  Likewise, cosine waves are necessary to accurately estimate the cosine components of your signal. Since most real signals are neither sin or cos waves, you need to use both to estimate them properly.

In addition to that, the DC component can only be represented by using a cosine wave (since sin(DC) = 0)

I personally don't think that you need to fully understand the math behind FT to realize that both waves are necessary.


If some values in $X[k]$ are zero, you cannot determine the impulse response of the complete channel bandwidth ($-f_\mathrm s / 2 \ldots f_\mathrm s / 2$, with sampling frequency $f_\mathrm s$). The IFFT of the estimated $\hat H[k]$ is rather a low pass filtered version of the channel. This might or might not be what you want.

In reality, you can never measure the complete bandwidth of a channel, as it is virtually infinitely large. For a communication system it is usually sufficient to know the channel transfer function of the band in which you intend to transmit a signal. If you transmit only on a fraction of subcarriers of an OFDM system you're effectively narrowing the transmission band and you might be satisfied with knowing the impulse response of this narrower band.


There was something like that in Aphex Twin – Windowlicker (1999 single), visible in the spectrogram with a logarithmic frequency axis:

According to the album's Wikipedia page, the image-to-audio conversion was made with Metasynth. I don't know the specifics of Metasynth, but a similar effect can be achieved by calculating the inverse discrete Fourier transform (IDFT) of each pixel column and concatenating the resulting time-domain frames. To make it sound smoother, the time-domain frames can be cross-faded and the phases of the frequency bins can be randomized to disperse the times when the phases of the frequencies match. The image can be distorted so that the vertical coordinate corresponds to a logarithmic frequency scale.

The problem could be stated as: How can I create an audio signal that produces a spectrogram that looks like the input image. To make a really good conversion, similar methods could be used as in time-scale/pitch modification. Some of those try to maintain phase coherence of frequencies (particularly of those that fall between frequency bins), across time-domain frames. Some of those try to maintain the time coherence of transients – things that are more concentrated in time domain than in frequency domain.


I beleive that h can be a vector or scalar .. depends on channel taps i.e. h

y[n] = hx[n] + w[n]  -- indicates  h convolved with x[n]  and then added to noise w[n]


A short burst of frequency f contains a ton of higher frequency spectra energy which is needed to make the burst short, or rectangular windowed, rather than infinitely long.  Thus, a finite burst no longer meets the Nyquust criteria for sampling at just a bit over 2 times f.

The transform of a rectangular window is a Sinc, which has infinite support in the frequency domain.  But a really really long window has a Sinc that decays faster, and thus may go below your noise floor soon enough.

If you want to sample at just a bit over 2 times f and not end up with aliasing, the signal may need to be stationary for a very very long time to allow reconstruction.

Or you can sample at a higher sample rate to acquire a sufficient amount of the spectra of the window that limits the length of a much shorter time domain burst.

So there's a trade-off between how close you can get to Nyquist frequencies versus the length of the signal.


Is it possible exctract sinusoids from non periodic signal?

Fourier's theorem say almost any (non-pathological) waveform, periodic or not, can be decomposed into sinusoids (or complex exponentials).  Whether, or how well, those sinusoids correspond to any underlying pseudo-periodic phenomena or not is another issue.

Note that with a DFT, you may need to interpolate a periodic signal between the DFT result bin sinusoids.


Is it possible the difference between OFDM-CP and OFDM-ZP about 10dB

I tried using MMSE equalizer, it worked well. single tap equalizer is not always working with ZP as well as CP OFDM.


How to find the normalized coefficients?

I'm not sure how you define your normalized difference equation, but in any case you can equate the $b_{ij}$ coefficents in the sos matrix with your $c$ coefficients, and the gain that tf2sos returns would then be your normalization constant $b$.

I think you already know, but just to be clear: every row of the sos matrix contains the coefficients for one second-order section, starting with the 3 numerator coefficients, and then the denominator coefficients (the first of which is always 1). You can incorporate the gain into the first section.


For Schwarz inequality for 2 signals $s_1(t)$ and $s_2(t)$, equality holds if $s_1(t) = c\cdot s_2(t)$ ($c$ is constant). Does $c$ have to be real?

Schwarz Inequality for continuous-time Complex valued functions is given as follows:
$$\left|\int^{\infty}_{-\infty}f(t)^* \cdot g(t) dt \right|^2 \le \int^{\infty}_{-\infty}\left|f(t)\right|^2dt \cdot \int^{\infty}_{-\infty}\left|g(t)\right|^2dt$$
As you can see that, on left hand side magnitude is taken after integration making it real number and the quantity on right hand side is also real number because magnitude-squared is taken before integrating. Thus, the comparison is valid and can be done for all Complex-valued continuous time functions which are Square integrable. Meaning the Inequality  exists for all functions $f(t)$ satisfying the following:
$$\int^{\infty}_{-\infty}\left|f(t)\right|^2dt = c < \infty$$


Calculate Max Number Of Harmonics for given Sample Rate?

No, you don't just multiply by 2, you multiply by every integer.  You can have a 2nd harmonic, a 3rd harmonic, 4th, etc.  (Only multiplying by 2s would tell you the number of octaves, not the number of harmonics.)

So take Nyquist frequency, and divide by the fundamental frequency

  22050 Hz / 20 Hz = 1102.5

Then round down to the nearest integer.  So you can fit 1102 harmonics.  The 1102th harmonic is 1102*20 = 22040 Hz, while the 1103th harmonic is 1103*20 = 22060 Hz, which is too high.

(Note that some consider the fundamental as the 0th harmonic while others consider it the first harmonic.  I am counting it as the first, 1*f.)


How to make semilogx to get a plot gain against log frequency with several values? [beginner]


Why shouldn't the magnitude of a practically realizable frequency selective filter be constant for a range of frequencies?

Note that all-pass filters with a perfectly constant magnitude response over the whole frequency range are realizable. A first-order discrete-time all-pass filter has the following transfer function:

$$H(z)=\frac{1-az^{-1}}{a^*-z^{-1}}\tag{1}$$

If $|a|>1$, the system $(1)$ is causal and stable. It is straightforward to show that

$$\left|H\left(e^{j\omega}\right)\right|=1,\qquad\forall\omega\tag{2}$$

i.e., the system's magnitude response is constant.

What is not possible is to have a constant magnitude response over only a part of the possible frequency range $[0,2\pi]$ (or, in continuous time, $(-\infty,\infty)$).

This can be seen as follows. First, design an ideal frequency-selective filter by computing the IDFT of some ideal frequency response (with a constant gain in the pass band(s), and a gain of zero in the stop band(s)). Let's consider discrete-time filters, but the discussion is basically the same for continuous-time filters. Let the resulting ideal impulse response be $h_{id}[n]$. First, we shift that impulse response to the right to have most of its energy at positive indices $n>0$:

$$\tilde{h}_{id}[n]=h_{id}[n-n_0],\qquad n_0>0\tag{3}$$

In order to obtain a causal filter, we need to multiply that shifted impulse response by a sequence $f[n]$ satisfying

$$f[n]=\begin{cases}0,&n<0\\1,&n>n_1\ge 0\end{cases}\tag{4}$$

with some non-negative $n_1$. The values of $f[n]$ between $n=0$ and $n=n_1$ can be chosen arbitrarily, but in practice we would choose them such that $f[n]$ increases monotonically in the interval $[0,n_1]$. The final impulse response of our causal filter is

$$h[n]=\tilde{h}_{id}[n]f[n]=h_{id}[n-n_0]f[n]\tag{5}$$

Note that $h[n]$ is still infinitely long.

In the frequency domain this corresponds to

$$H(e^{j\omega})=\left(H_{id}(e^{j\omega})e^{-jn_0\omega}\right)\star F(e^{j\omega})\tag{6}$$

where $\star$ denotes convolution. This convolution destroys the piecewise constant property of the original frequency response $H_{id}(e^{j\omega})$, and the magnitude of the resulting causal frequency response $H(e^{j\omega})$ is neither piecewise constant, nor does it have infinitely sharp transitions.


Detect marker orientation inside an image

In essence, you are looking for a way to do planar object detection.

Well, the most basic approach is to find a correspondence between your planar template and the image using a feature extractor, such as SIFT. You could then cast the detection problem as finding the subset of points, which have similar descriptors. OpenCV has that implemented already:
http://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_feature2d/py_feature_homography/py_feature_homography.html

http://docs.opencv.org/trunk/dc/dc3/tutorial_py_matcher.html

The next thing is to use template matching kind of approaches. For severe viewpoint this might be hard to achieve, but there are ways, such as this one.

If you can also design a marker in the image space, then there are approaches which let you to read a randomly generated marker. For example:

  T. Birdal, I. Dobryden, S. Ilic  X-Tag: A Fiducial Tag for Flexible
  and Accurate Bundle Adjustment  International Conference on 3DVision
  (3DV), Stanford University, California, USA, October 2016
  http://campar.in.tum.de/pub/tbirdal20163dv/tbirdal20163dv.pdf

Take a look at the related work section of the references I gave here. This could broaden your knowledge in the field. You might also like to read a survey in the field, such as:
http://www.irisa.fr/lagadic/pdf/2012_fcv_uchiyama.pdf


Without knowing what you want to get, I will assume you are interested in periodicities of the elevated dot patterns. They can give you access to roughness measures, texture attributes.

Visually, 2D periodicities are apparent from your data point $(x,y)$ locations. At different scales though: beside a potential periodic base pattern, little streaks of 5 to 10 dots, closely spaced, can be observed.

The situation will depend on how to want to treat the data. If you treat it as an image, putting white values (or whatever) between dots, you are already doing interpolation. In that case, a true 2D DFT would take a separable 1D DFT on both rows AND columns. So, you definitely need to venture into 2D, unless you will only  see horizontal frequencies; look at the last line of your image: there are 8 or 9 periods only, you are missing out the small slope.

I   cannot propose data interpolation without knowledge of the underlying process, and the $z$ value. My take is that you ought to try with true 2D non-uniform DFT first. Interpolation causes smearing.

In recent question (what is $1/f$ property of the visual world?), Frequency Synthesis of Landscapes was evoked, and a parametric spectrum modeling (eg decay in $1/f^p$) can provide you with parameters to compare different surfaces. For some spectral shapes, this can be naturally done in a Fourier domain. There should be related techniques under the field of "Surface roughness analysis".

My basic thinking for a procedure is the following: your different scans of surfaces may have different $(x,y)$ locations, and sizes. So it is not easy to compare their properties directly in the space domain, unless you regrid, interpolate, etc. which can cause harm. Non regular Fourier transforms can regularize the outcomes: from any two scattered grids, you can cast Fourier coefficients to the same 2D regular frequency grid, and infer global parameters from them: peaks, moments, parametric fits, etc.

So you can check NUFFT (NFFT, USFFT) Software and the related keywords:

  algorithms have been developed to overcome this limitation - generally
  referred to as non-uniform FFTs (NUFFT), non-equispaced FFTs (NFFT) or
  unequally-spaced FFTs (USFFT)

with a related paper: Accelerating the Nonuniform Fast Fourier Transform, by Leslie Greengard and June-Yub Lee. Apparently there are Matlab versions too:
NUFFT, NFFT, USFFT.

If this is not enough, you can try (lifting) wavelet transforms on 2D irregular grids, to catch the proper scale and local details. You could mesh it with MeshLab. And then use lifting wavelets for a multiresolution vision.


Is the book wrong? -- How to sketch frequency response obtained from H(z)?

As far as I can tell, the book is not wrong. (By the way, you're studying a 3-point moving average filter here, and it's useful to know about such filters.)

Your 1st question: That center dot on the Mag. and Phase curves serves no purpose that I can see. I suggest you not worry about that center dot.

Your 2nd question: The Phase curve is generated by evaluating the angle of the  $e^{-j\hat{\omega}}$ phase term, in the last line of the $H(\hat{\omega})$ expression in Part (c) of the Solution, as $\hat{\omega}$ goes from $-\pi$ to +$\pi$. As $\hat{\omega}$ goes from $-\pi$ to +$\pi$, the Solution's $\hat{\omega}/2\pi$ Phase curve freq axis label goes from -0.5 to +0.5.

(Note: as you study the Solution know that the $H(z)$ in Part (c) is not the $H(z)$ in Part (b). The $H(z)$ in Part (c) is the $H(z)$ from the Solution's Part (a).)

Your 3rd question: It's a good question. In the literature of DSP there are a number of different notations (different symbols with different dimensions) used to represent frequency. And those various notations certainly can be confusing at first. mmswe, have a look at the following web page, it might help you with your 3rd question:

http://www.dsprelated.com/showarticle/199.php


How to show that the signal $x_n = A\cos(\omega n)$ can be fully predicted by a system with two weights $w_1,w_2$

HINT:

Use

$$2\cos(x)\cos(y)=\cos(x+y)+\cos(x-y)$$

with $x=(n-1)\omega$ and $y=\omega$.


Error implementing Harmonic Product Spectrum algorithm

Pad the high frequency portions of all the down-sampled vectors with either zeros, or some noise floor estimate, depending on your desired frequency search range, to make them all the same length or dimension.


Why does decimation make a system time variant?

HINT:
Think about the following situation: you have an input signal with alternating ones and zeros, and you decimate by a factor of $2$. Now shift the input signal by one sample and compute a new output signal by decimating by a factor of $2$. Is the output signal a shifted version of the previous output signal? If the system is time-invariant, shifting the input signal causes the output signal to shift by the same amount.


Fourier Spectra : Significance of the Negative Amplitude

Note that in general the Fourier transform of a function is a complex-valued function, so in general it is not only positive or negative.

Roughly speaking, the magnitude of the Fourier transform says something about the presence of certain frequencies components in a signal, regardless of the phase (or sign, in the real-valued case). The phase determines the temporal alignment of the different frequency components. So a real-valued negative value of the Fourier transform just implies a 180 degrees phase shift of the respective components.


when I read the original autotune patent few steps made me understand that everything was done in time domain (in the past, I don't know today), they didn't mention anything about overlap and add, it made me wonder if the pitch detector is so good that they didn't have to overlap and add, could they always skip or add periods in the exact position? (just out of curiosity I ended up implementing the pitch detector described in this patent, and honestly that's not all I imagined lol)

Take a look in the red marked piece of the patent paper:

Of course, the patent description also says nothing about keeping the formants, but to me it is clear that the above description is just a time domain method to apply time scale and change the pitch using resample.

By the way, I wrote a pure time domain "pitch correction" keeping the formants based in the papers from @robertbristow-johnson and Keith Lent

youtube demo from my code here:

https://www.youtube.com/watch?v=6ns5K1FHtd4


What is meant by "spectral whitening" in DSP?

Spectral whitening is usually an attempt to make the spectrum of the signal "more uniform".  One reason this might be a good thing to do is that it can have the effect of making the autocorrelation of the signal "narrower" (and closer to a Dirac delta, for discrete-time signals). This can help localize in time.

  What effect does spectral whitening have when used in image processing? (visually or otherwise...)

It's generally not pretty.  Most images are "low pass" (most of the information is in the low frequency part of the spectrum).  One simplistic approach to whitening in images is to do a column-wise (or row-wise) difference (i.e. diff in matlab).

This will mean negative pixel values, which generally do not map to anything sensible with standard images.

This example shows how prewhitening can improve localization in image processing template matching. The picture from that link is:

  Where might spectral whitening be useful in audio processing or analysis?

If you are trying to localize (in time) the onset of a sound, then it's possible that spectral whitening can improve this.  It's also possible that it can reduce (disimprove) the SNR.

  What would a spectrally whitened audio signal sound like?

For audio of speech or music, it will tend to bring in more higher frequencies.


Why Butterworth filter always starts my signal from from zero mark (amplitude)?

Filters do have a delay (a lag) since they do not act immediately on your signal. Also all samples before the time 0 are zeros, thus in general you will start from the "zero mark", as you said (just imagine your filter equation with all zeros).

There are ways to make a filter to have a zero lag. It is done by so called zero-phase filtering, also known as forward-backward filtering. The way you do it is by filtering your signal twice - in forward and in backward direction. Obviously this can work only for offline applications. In MATLAB you can do it very easily using the filtfilt function.

Please find below a code corresponding to your case.

fs = 100;
T = 2;
t = 0:1/fs:T;
f = 2;
s = 1.5*cos(2*pi*f*t) + 0.7*sin(2*pi*33*t);

[B,A]=butter(6, 5/fs,'low');
S=filter(B, A, s);
Sff=filtfilt(B, A, s);

plot(t, s)
hold on
plot(t, S)
plot(t, Sff)
grid on
legend({'Original signal', 'Signal filtered using filter', 'Signal filtered using filtfilt'})

That yields what I suspect you want:


How do I make sense of the cosine wave having Fourier Transform coefficients which have infinite magnitude?

Your function is both integrable, and square integrable, because it is continuous on the finite support $[0\;2\pi]$. In this context ($L_1 \cap L_2$), its Fourier transform is pretty well defined as a classical function (not a distribution like the Dirac), and continuous. So, your formula is incorrect.

So there is no infinite magnitude. Moreover, one cannot simply say that Dirac have "infinite magnitude".

Yet, this would not imply that a perfect cosine wave can not be produced in real life; produced by whom? However, there is little chance that somebody could observe, or capture, a perfect cosine wave in real life (caveat: if real life can be described by a continuous time variable, and amplitude could be continuous as well).


The Gabor wavelet is (almost) the same as the Morlet wavelet (some authors distinguish these two by an additional constant to fulfill the wavelet's admissibility condition). And yes, it is the same function as in the Gabor transform (gaussian windowed oscillation).

The application of phase coherence should be independent of the underlying wavelet and the nature of the transform (Wavelet, Gabor/STFT, other weird metrics). Essentially important is the concept of reproducing kernels and the intrinsic (Heisenberg) uncertainty that influence you measure of phase coherence and must be taken care of.


Q1: Dependence of estimation performance on number of data points

Since LMS and RLS are adaptive filters, their estimation performance improves as the number of their iterations increase. Hence more data points will make their outputs closer to the expected performance, until the convergence is achieved (this requires either a WSS data or a slowly changing, statistically,  nonstationary data, so that filter can track it's statistical character). Once the filter is operating in convergence conditions, then there won't be any improvements in its output, providing what is called as the (minimum) steady-state error. This steady-state error depends on a number of factors but not on the data length.

But until the convergence is achieved, the estimation error decreases as more data points processed. However most typical adaptive filter users would be interested in its steady state (convergent) results rather than the transient response.

Note that increased data points can either come from a longer observation or from higher sampling rate, the results of which can be different. For example for a Kalman filter (in its extended mode) with mechanical applications, an increase in sampling rate can actually improve the steady state error as well.

  Q2: Dependence of convergence on number of data points

For both LMS, RLS and Kalman filters, convergence primarily depends on the number of data points being processed; i.e., number of iterations. However this convergence rate is different for different filter types, reflecting their complexity and/or sophistication. The simple LMS filter has the slowest rate of convergence (roughly 10 to 20 times the filter tap weight length) whereas the RLS and Kalman filters display a convergence rate of roughly $2$ times the filter length (Haykin_Adaptive Filter Theory) for WSS inputs.


Ultrasonic speaker calibration and emitting calibrated signals

Found the answer after looking at the references that Jim Clay mentioned in the comments, thank you Jim.

I made the mistake of only considering the magnitude which results in zero-phased signal and can not sensibly be used for emission, at least not in this setup.

The code I finally ended up using can be seen below.

The script adheres to the naming convention of keeping time domain signals in lower-case and frequency domain signals in upper-case.

% Align and sum all files called Mandag*
files = dir('Mandag*');

% Where in the recordings the signal is
rng = [1.5e6, 1.52e6];

% Initialize the xh vector
[xh, fs] = wavread(files(1).name, rng);
xh = xh(:,1);

for i=2:length(files)
    y = wavread(files(i).name, rng);
    y = y(:,1);
    % Determine offset between xh and y
    [~, off] = max(xcorr(xh', y'));
    off = length(xh) - off;
    % Shift signal appropriately
    if(off < 0)
        y = [zeros(1, -off), y(1:end+off)'];
    elseif (off > 0)
        y = [y(off:end)', zeros(1, off-1)];
    end
    xh = xh + y';
end

% Average
xh = xh/length(files);

% Location of the 20ms signal
xh = xh(2306:12306-1);

% Normalize
xh = xh / max(xh);

% Apply a moving average filter on xh to reduce noise. Window size of 4 was
% experimentally determined to give the best results
n = 4;
B = zeros(n, 1);
for i=1:n
  B(i) = 1/n;
end
xh = filter(B, 1, xh);
xh = xh / max(xh);

x = wavread('sweep.wav');
x = x(1:2:end);            % Sweep generated @ 1MHz, decimate
                           % to have same length as xh

% Transform x into frequency domain and determine H
X = fft(x);
H = fft(xh) ./ X;

% Vector indices to choose only frequencies of interest
starti =  20e3 / 50;
endi   = 100e3 / 50;
rng    = starti:endi;
irng   = (length(x) - endi) : (length(x) - starti);

% Zero out unwanted frequencies
X = [zeros(1,      starti - 1   ), X( rng)', zeros(1, length(X)/2 - endi) ...
     zeros(1, length(X)/2 - endi), X(irng)', zeros(1,      starti - 1   )]';

% Deconvolve x with h
X_deconv_H = X ./ H;

% Transform X/H to time domain and normalise
x_deconv_h = real(ifft(X_deconv_H));
x_deconv_h = x_deconv_h / max(x_deconv_h);

% Save the deconvolved sweep
wavwrite(x_deconv_h, fs, 'deconvolved_sweep.wav');

% Generate  spectrograms of xh and x_deconv_h
winsize = 512;
overlap = round(.99 * winsize);
figure(1)
specgram(xh, winsize, fs, hann(winsize), overlap)
colorbar
figure(2)
specgram(x_deconv_h, winsize, fs, hann(winsize), overlap)
colorbar

The spectrograms of x conv h and x deconv h can be seen below:

These seem plausible to me although there is some noise in the deconvolved signal.

Next test will be to see if emitting x_deconv_y gives something resembling x barring those frequencies that the speaker cannot emit.

Update with test results

We redid the measurements described above using a logarithmic down-sweep. These results seem to suggest that the method works.

The verification test consisted of emitting X / H and expecting to get X back, i.e. equal energy in all frequencies.
As the worst output frequency is about 20dB weaker than the best, the highest output level is expected to be that much lower.

The signal that was emitted:

The time series and spectrogram of the recorded signal look like this:


For linear filters, binomial filters, which can be considered as FIR approximation to Gaussians, are very simple. Their coefficients are given by:
$$ h_k=\frac{ \binom{n}{k}}{2^n} \,.$$
For instance, you can get $H_3 = [ \frac{1}{4},\frac{1}{2},\frac{1}{4}]$, or $H_5 = [ \frac{1}{16},\frac{1}{4},\frac{3}{8},\frac{1}{4},\frac{1}{16}]$. Since their coefficients are positive, they do not overshoot.

Here are graphs for $n = 5,17,25,31$:

Savitzky-Golay filters, performing polynomial interpolation, are interesting too, but tend to   overshoot with higher degrees.

Do not fool yourself though: as long as the filters are linear, there is a notion of frequency very close behind.


Does audio upsampling create noise/artifacts or degrade the signal?

The up-sampling process will always change the signal in some measurable way. However, if it's done properly the changes are negligible and don't result it any audible difference. Most commercially sample rate converters (hardware or software implementations), do a really good job at this.

Off course, if done badly, upsampling can result in clearly audible signal degradation. I'm not familiar with Apple's implementation but I would assume that they got this correctly.


From that same page you link to:

  A signal is said to be oversampled by a factor of N if it is sampled at N times the Nyquist rate.

So, you're sampling at $1\cdot10^4\,\text{Hz}$, that's $N=\frac{10^4}{140}\approx 71$; no need to subtract anything.

  Does this oversampling only reduce the noise introduced from the output of the analogue sensor (given the low bandwidth of the sensor) to the end of the ADC process?

Oversampling alone doesn't reduce noise. It's what you do afterwards with the signal. Oversampling can be necessary to e.g. to be able to filter out the noise that can be shaped to your advantage.

  How much can I reduce the oversampling (which introduces a processor overhead) by for a negligible increase in the amount of noise to the signal?

You need to define "negligible"; also, you need to think about the noise is distributed in frequency domain during your sampling process; there's no general answer.

  Could I do anything else to reduce the amount of oversampling, e.g increase the ADC resolution?

I'm confused: The only way to reduce oversampling is to reduce the sampling rate.

Maybe you mean

  to reduce the necessary amount of oversampling

Well, that fully depends on your signal, and what you want to do with it, and you forgot to tell us about either!


Essentially, OS is slightly more efficient since it does not require the addition of the overlapping transients.  However, you may want to use OA if you need to reuse the FFTs with zero-padding rather than repeated samples.

Here is a quick overview from an article I wrote a while ago

  Fast convolution refers to the blockwise use of circular convolution
  to accomplish linear convolution. Fast convolution can be accomplished
  by OA or OS methods. OS is also known as “overlap- scrap” . In OA
  filtering, each signal data block contains only as many samples as
  allows circular convolution to be equivalent to linear convolution.
  The signal data block is zero-padded prior to the FFT to prevent the
  filter impulse response from “wrapping around” the end of the
  sequence. OA filtering adds the input-on transient from one block with
  the input-off transient from the previous block. In OS filtering,
  shown in Figure 1, no zero-padding is performed on the input data,
  thus the circular convolution is not equivalent to linear convolution.
  The portions that “wrap around” are useless and discarded. To
  compensate for this, the last part of the previous input block is used
  as the beginning of the next block. OS requires no addition of
  transients, making it faster than OA.


Extracting the DC component of a signal, that proportional to the displacement about the y-axis

DC component is the zero frequency component. If you look at the DFT expression, this would mean summing up the input. So in mathematical terms this would be the sum (or the running sum) based on your application.

In terms of filtering again this would be an averaging filter. Like moving average.


Why edge sharpening produces high frequency?

Edges are not the best defined features in images.

However, they can be associated, locally, at a certain scale, with relative variations in intensity along a first direction, combined with a relative smoothness in a complementary (for instance orthogonal) second direction. When looking along the first direction, the 1D intensity profile exhibit variations, faster  than in the other direction. Several typical edges profiles exist, as depicted below.

Those are typically enhanced by derivative-like operators, which are weakly sensitive to flat, linear or generally regular intensity variations. In other words, they behave as band-pass or high-pass operators.
 In the other direction, gradient detectors offer more smoothing (more low-pass).

Another vision is to smooth an image, and subtract this result  from the original image. By removing low-frequencies, higher ones are relatively amplified. But to avoid noise amplification, some smoothing is often included in the process.

Hence, linear or not, edge operators tend to emphasize mid- to high frequencies relatively to the lower ones. However, all this may happen at different scales, depending on the edgy objects of the picture.

Now,  Is this new information reliable?

Not really, if you rely on mere filtering or enhancement. It is easy to produce spurious edges. But within a restoration or inverse problem context, with models and an objective function to minimize, more reliable information can be obtained, as long as the models mimic well the physics of your acquisition scheme.


The IEEE signal processing literature deviates somewhat from the terminology used in statistics and actually from itself so your confusion is justified.
If $\mathbf{x_k}$ is a vector signal the quantity at sample time $\mathbf{k}$,

$$ \mathbf{R}=
\sum_k \mathbf{x_k \, x_k^H}
$$

is often called either.

The SP transactions historically focused on application areas like like radar, sonar, speech, audio, telephony, .... where the signal propagates in a way where dc isn’t propagated.  The long term time average of amplitude tends to zero.

Music (assumptions) works on 2nd order statistics and is stationary over some interval so implicitly the signal is random.  Under these circumstances the average values of the data will be zero. The algorithm is generally used in array processing.

Typically in IEEE SP papers covariance denotes that the mean vector is subtracted out so the correlation (as shown above) of a zero mean signal vector and the covariance is the same thing.

The interesting thing about algorithms that use 2nd order statistics is that they often work well on deterministic signals.

Authors in Kalman Filtering papers are more specific.

In Statistics, the Pearson correlation coefficient is not the same as the usage in array processing.

There are more than a few words that have different but similar meaning, like coherence.

Even in SP, the terms mean, average, and expected value are often used in ways that lead to confusion.


Signal in frequency domain with OpenCV dft

Your sine wave should use 2pi instead of frame rate.  Also, your frequency data should have window size/2 points and go from 0 to frame rate/2, not frame rate.


Finding the $\mathcal Z$-transform of $((-\frac{1}{3})^n + \frac{1}{2})^n \mu[n-2]$?

Since $\left| -\frac{1}{3}\right|< \frac{1}{2}$ and $n$ starts at $2$, one can think about an approximation based on the first terms of the polynomial development of $a_n = \left(\frac{1}{2}\right)^n\left(1+2\left(-\frac{1}{3}\right)^n\right)^n $. Using a binomial (polynomial) development up to degree $2$ on the second term only, one gets:

$$a_n \approx b_n= \left(\frac{1}{2}\right)^n + 2n\left(-\frac{1}{6}\right)^n + 2n(n-1)\left(\frac{1}{18}\right)^n \,.$$

Details: since coefficients are zero for $n<2$, for $n=2$,

$$a_2=\left(1+2\left(-\frac{1}{3}\right)^2\right)^2=  b_2 = 1+ {2 \choose 1}\left(2\left(-\frac{1}{3}\right)\right)+ {2 \choose 2}\left(2\left(-\frac{1}{3}\right)\right)^2$$

exactly. Then for $n>2$, the term $c_n=2\left(-\frac{1}{3}\right)^n$ tends to be small. You can approximate $(1+c_n)^n$ by  $1+ {n \choose 1}c_n+ {n \choose 2}c_n^2$.

The following graphs provide an evaluation of the approximation, up to $n=8$. Top graph superimposes $a_n$,  @MaximilianMatthé's approximation (MM), and what we get using the first (D0), the first two (D1), or the three terms (D2) of the above series. $D0$ (only $\left(\frac{1}{2}\right)^n$) is clearly insufficient. The others are quite close. The bottom graph shows the absolute approximation errors. As you can see, $D1$ and $D2$ converge very fast when $n$ grows, and are already quite good for the first terms.

The nice thing with these approximations is that you not only reuse the classical properties of  the $\mathcal{Z}$-transform: linearity, time-shift and scaling, but also differentiation (up to order $2$), because of the factors $n$ and $n(n-1)$:

$$ n x_n \to -z\frac{dX(z)}{dz}\,.$$

A Matlab code for reproduction:

n = (2:8)';
s = ((-1/3).^n+1/2).^n;
s0 = ((-1/3).^n+(1/2).^n);
sa = ((1/2).^n);
sb = ((1/2).^n).*(1+2.*n.*(-1/3).^n);
sc = (1/2).^n+2.*n.*(-1/6).^n + 2.*n.*(n-1).*((1/18).^n);

subplot(2,1,1)
plot(n,[s,s0,sa,sb,sc],'.-');axis([0 max(n) -Inf Inf]); grid on
legend('Orig.','MM','D0','D1','D2')
h=xlabel('$n$ index');set(h,'INterpreter','latex');title('Series')
subplot(2,1,2)
semilogy(n,abs([nan(size(s)),s0-s,sa-s,sb-s,sc-s]),'.-');axis([0 max(n) -Inf Inf]); grid on
legend('Orig.','MM','D0','D1','D2')
h=xlabel('$n$ index');set(h,'INterpreter','latex');title('Absolute approximation errors')


When do we divide by $2\pi$ while solving Nyquist problems?

Well!  all signals in this world are made up of sum of different rotations(sinusoidals) - different in three senses:
a. how big is the amplitude (A)
b. how fast is the rotation ($\omega$)
c. where is the starting point of the rotation (phase $\phi$)
Fourier made this very clear.
How do we measure rapidness of the rotating signals(sinusoidals) : by their angular velocity $\omega$ which is given in radians/seconds. This is the correct and most appropriate unit for measuring rotations - how much of angle is being covered per unit time! If a constituent signal is rotating by $200\pi$ radians in 1 second, we say $\omega = 200\pi$ rad/sec.
But we can also measure the rapidness of these rotating signals as number of rotations per unit time. That is $f$ expressed in $Hertz$ : number of rotations per second.
How are $f$ and $\omega$ related? One complete rotation covers a full circle meaning an angle of $360^o = 2\pi \ radians$. That means if the signal is making $200$ rotations per second then it is covering an angle of $200*2\pi = 400\pi$ per second. So, the relation between $f$ and $\omega$ is basically:
$$\omega = 2\pi f$$
That is why you divide $\omega$ by $2\pi$ when you want to express the frequency in $Hz$.
The example you have given is $\cos{\frac{\pi t}{3}} = \cos{\frac{2\pi t}{6}}$. The period of this sinusoidal is $T = 6$sec, therefore, the frequency in $Hz$ will be $f = \frac{1}{6}$ and in radians/sec will be $\omega = 2\pi f = \frac{2\pi}{6}$ rad/sec.
(Think about why the period of this sinusoidal is $6sec$. Figure out that the sinusoidal will repeat after $t = 6sec$. Figure out that the sinusoidal makes 1 full rotation of angle $2\pi$radians in $6sec$.)
And the sampling frequency is twice the frequency of the sinusoidal giving :
$$f_{sampling} = \frac{2}{6} = \frac{1}{3}sec$$
Meaning Sampling period is $\frac{1}{f_{sampling}} = 3sec$.