Morphological operations to isolate and preserve network graph structure? Suppose there is an image with long, connected lines which form connected, closed cycles in an image.
Is there a standard morphological approach to isolate only the lines, and drop all other pixels?

I can imagine a couple of slow approaches:

convolutional approach where a window -- say, 50x50 -- is passed along the image, and any non-zero pixels which are not connected to a pixel on the perimeter of the window are set to 0.

Skeletonize - threshold the components - close - skeletonize ...: the skeleton is calculated, small components are dropped, disconnected skeletons are closed, and the skeleton is recalculated and re-thresholded. Repeat.

But it seems like this would be a common problem, and there should be an easy solution...?

OpenCV offers a set of routines embracing all the operations necessary for graph network extraction from pixel images. Yet, to my knowledge, there is no a "graph network extraction" toolbox in OpenCV on the manner of MATLAB toolboxes that supply integrated solutions for DSP, image processing, etc. tasks.
The search for a solution of the kind brought me a Nature article NEFI: Network Extraction From Images by M. Dirnberger, T. Kehl & A. Neumann.
The NEFI tool is freely available, both in the form of release binaries and the source code. I have no practical experience with this tool and cannot recommend it to you, but I like the Nature article and follow their guidelines to summarize the review of graph network extraction techniques. In fact, I've compiled my answer from the article text, attaching comments and references for your convenience.

A typical pipeline of a solution for graph network extraction combines algorithms from up to four different classes: preprocessing, segmentation, graph detection and graph filtering.
PREPROCESSING
Uses operations on pixel images, like filtering (blurring, denoising) and maybe color adjustment, to ease the operation and, to a degree, to improve the results of a segmentation step.
SEGMENTATION
Separates the image foreground, i.e. the structures of interest representing a graph network, from the remaining image.
It seems that the first (historically) segmentation method is Otsu's method (adaptive thresholding): a threshold that gives the best separation of "background" and "graph object" pixel classes in gray levels would be the best threshold.
Watershed algorithms are based on the "flooding" process performed on the gradient image (priority-flood). The algorithms define and use a set of markers, pixels where the flooding shall start. Implementation and application: see the article On Graph Extraction from Image data.
Otsu's and watershed segmentation approaches work on grayscale images.
The GrabCut algorithm based on graph cuts works on color images and, similar to watershed algos, also receive an additional input from markers. The algorithm is implemented in an openCV function
cv.grabCut (image, mask, rect, bgdModel, fgdModel, iterCount, mode = cv.GC_EVAL)

, see an OpenCV tutorial Foreground Extraction using GrabCut.
See also an CNN implementation: https://grlearning.github.io/papers/23.pdf
GRAPH DETECTION followed by GRAPH FILTERING
Morphological operations (thinning) performed on segmentation results give a skeleton of the segmented image. The skeleton is used to find the graph node positions and consequently the graph edges.
If you've read up to this point, you may agree that a one-piece morphological operation of your description,

a standard morphological approach to isolate only the lines, and drop
all other pixels

, has yet to be developed.
To the moment, in order "to isolate only the lines, and drop all other pixels", the software developers can only offer a pipeline of at least two stages: segmentation (e.g., with a $cv.grabCut(..)$ function) and graph detection (thinning procedure).
The output of graph detection algorithms is further improved in the graph filtering stage using heuristic or user assisted methods. The NEFI creators developed an approach exploiting the structure of the extracted graph.
The graph filter operation can also be considered partially overlapping with your item 2 description.

I hope this answer would help you develop your own solution, or you will decide to contribute to the NEFI project or the other graph extraction projects on GitHub.


How to do convolution in Fourier Series? Two signals are given to me :

$$x(t)=\cos(4\pi t)$$
$$y(t)=\sin(4\pi t)$$

I have founded their coefficients as follows:

$$a_k = a_1=a_{-1}=\frac{1}{2}           $$

$$b_k = b_1=b^*_{-1}=\frac{1}{2j}           $$

Now I am stuck at convolution. In the formula of the convolution:

$$\sum_{l=-\infty}^{\infty}a_lb_{k-l} $$

What are these $l$ and $k$ and from the coefficients I found what is their respective value?

Following your notation, if we define the coefficients of the Fourier series of $z(t)=x(t)y(t)$ as $c_k$:

$$c_k=\sum_{l=-\infty}^{\infty}a_lb_{k-l}$$

Notice that the $k$ in which you evaluate the left side of the equality has to be the same as the $k$ in the summation. Therefore, for example, if you want to find $c_0$, you should calculate:

$$c_0=\sum_{l=-\infty}^{\infty}a_lb_{-l}$$

Then you have to calculate that summation. Given that $a_l$ and $b_{-l}$ are non-zero only if $l$ equals $1$ or $-1$:

$$\begin{align}
c_0
&=\sum_{l=-\infty}^{\infty}a_lb_{-l} \\
&=a_{-1}b_1 + a_1b_{-1} \\
&=\frac12 \frac{1}{2j} - \frac12 \frac{1}{2j}=0
\end{align}$$

I think this might have clarified the notation for you. You can go on and calculate the rest of the coefficients $c_k$ for all values of $k$.


How to visualize high-pass filtered images? When I try visualizing a high-pass filtered image, all I see is gray, similar to the middle subfigure in the attached figure from a related paper. The authors claim they normalize the image to have unit variance, so I tried plotting $(x - \mu) / \sigma$, which still resulted in a gray image. Here is my process for high-pass filtering:

Shifted DFT
Apply radial mask, zero out low frequencies
Inverse shifted DFT
Normalize image to $[0, 1]$

First of all beware of the difference between image sharpening and highpass filtering. The former is actually a high frequency amplified image, while the latter removes the low frequency and completely eliminates the DC component.

So assuming that you indeed wanted to remove low frequency content by appling a highpass filter, then your output will be a ghost like image with black background that shows image edges.

Highpass filtering can produce out of range values, compared to the original valid range. This can be corrected with the following:

Assuming valid data range as [0,1], then the following assures that it will remain in [0,1] again:

x[m,n] = (x[m,n] - x_min)/(x_max-xmin)

Assuming valid data range as [0,255], then the following assures that it will remain in [0,255] again:

x[m,n] = 255*(x[m,n] - x_min)/(x_max-xmin)

More generally to map a data range from [a,b], into [c,d], the following mapping can be used:

x[m,n] = (d-c)*(x[m,n] - a)/(b-a) + c

Note that a range correction may not always yield a best looking image. You might as well need to adjust distribution of intensities in nonlinear ways to pop out requested features...


Basic method for 2x oversampling? I would like to know if the following method of 2x oversampling is correct:

Interpolate:

Take an original signal sampled at 44100Hz as input
Upsample by adding a zero after each original sample to get a signal twice as long
Filter the new signal using a low pass filter (cutoff of 44100Hz based on new sample rate of 88200Hz)

Use new signal:

Do some additional filtering of the the new signal at 88200Hz

Decimate:

Take the final signal at 88200Hz as input
Filter using a low pass filter (cutoff of 44100Hz based on sample rate of 88200Hz)
Downsample by removing every other (odd) sample
Now have a signal that is 44100Hz again

What I would like to check is if the method is correct and my cut off frequencies for my low pass filters are right.

Interpolation (sampling frequency 44.1 kHz ➔ 88.2 kHz)

Your original 44.1 kHz sampled signal has frequencies up to 22.05 kHz, so you should lowpass filter at 22.05 kHz after dilution with zeros. Your filter should have a gain of 2. Otherwise the signal amplitude drops to half because you set half of the samples to zero. Like Jim Clay says, you can combine these things; there will be plenty of opportunity for optimization.

Downsampling (sampling frequency 88.2 kHz ➔ 44.1 kHz)

If your 88.2 kHz sampled signal is the one you lowpass filtered at 22.05 kHz, you can simply decimate by throwing away every second sample. Otherwise, first lowpass filter at 22.05 kHz to prevent aliasing, and then decimate. If you combine lowpass filtering and decimation, you can optimize again.


Ambiguity function of multiple sines I was reading through the Ambiguity function, and I was wondering if it's possible to calculate analytically (or numerically) the ambiguity function of a sum of sine waves. In case of a single sine wave the envelope is a rect(), hence the ambiguity function is a $sinc()$ with variables the delay time and the Doppler frequency. Here's  the respective 3D plot. Is it possible to perform a similar reasoning for multiple sines?

Sum of sines isn't a concrete waveform, there is no general form of its ambiguity function. But for any given sum of sines you can calculate this function by evaluating cross correlation for various time and Doppler deviations. Though one rule exists, if the number of sinewaves in composition grows, the ambiguity function will become narrower in both frequency and time dimensions, since such a signal will be more and more noise like.

Since autocorrelation isn't linear function of some signal $s$ (its quadratic as you know cause it is sum of multiplications of $s$ and its shifted copy), so you can't express autocorrelation and more generally ambiguity function of sines composition as a sum of their autocorrelations (or amb. func-s.). Actually it is the sum of entering sines' autocorrelation functions and all possible entering sines' cross correlation functions. It could be something huge. So the most general way is to compute of such a function for given sines sum directly with knowing formula of autocorrelation.


Automatic Gain Control for voice I am tasked with designing the logic for an automatic gain control. It won't be a traditional AGC in that it has to run on a computer and continuously modifying the volume is not possible in the same way as it is in hardware. It will be more of a "roughly set the input gain for me automatically if it's way off" than an AGC per se. However, I plan to approach the design in the same way.

This will be used primarily for voice/speech, so what I want to do is prevent incoming silence (when the speaker is not speaking) from causing the gain to be set too high. I would like to avoid anything as complex as voice activity detection as well as complex ad-hoc methods that involve deriving more than a few simple parameters -- at least for now.

Are there any simple tricks I can use to to help with this problem, or is the best approach to use full-fledged voice activity detection?

Two things come to mind.

First, let's address the low-amplitude noise problem. This one is a typical textbook case for a noise gate. To give you an extreme example, these are sometimes even used on guitar distortion effects which are essentially insane amounts of controlled overdrive. I have to admit that I failed to find any literature on creating digital noise gates, but they exist abundantly as plugins so maybe digging through some open-source projects will score you a good basis for learning how to implement them. In essence, you measure windowed energy of a signal, and if it falls below a certain threshold, you fade your gain to zero (usually decaying exponentials are used). In any case, I'm sure you know what a noise gate is and I'm beating around the bush here. Make one and use it.

Secondly, when you say roughly set the input gain for me automatically if it's way off, that sounds like a textbook case for using a fuzzy controller. I'm linking to the Wikipedia page because I'm feeling a bit lazy, but you can find tons of great papers out there just on the first ten pages of Google results. The concept is exactly as you described it. You'll have several (say five for the sake of a simple example) levels of how close the gain is to what you think it should be:

- Way too low
- Low
- About right
- High
- Way too high

Based on these you will have some sort of fuzzy response like:

- Increase at a high rate
- Increase at a low rate
- Do nothing
- Decrease at a low rate
- Decrease at a high rate

These are my thoughts on how I would do it. You still need some means of measuring your fuzzy control input (assessing how good current gain level is), but that's something for you to work out. I'd say that if you're riding very close to clipping, that's High, if you've clipped, that's Way too high, and so on. Depending on what the spec considers acceptable, fuzzy control can sometimes let you get away with very simple solutions.


Is there a transformation filter to decode light signal through Glass Bricks? A friend is doing renovations and getting his entire ground-floor street-side wall replaced with glass bricks.  I told him that with all those bricks giving similar distorted views of the same room, and scattering light in similar ways, there must be some statistical technique to inverse filter the light that comes through to clean the signal and see what's behind the bricks---and that he should watch out or at least consider curtains.

Was I right?  Or is this impossible (or NSA level stuff), or can this be done with standard image processing libraries and a little stats?

In theory this is possible, in practice probably very hard, but I wouldn't say impossible. See a recent article by the Silberberg group for example, that demonstrates imaging through 'thin' turbid layers using scattered incoherent light. Also, search for works pioneered by A. P. Mosk on that matter (for example this one and this one). The main challenges are algorithm speed, or, by the time you'll find the correct transmission matrix, it'll change due to thermal fluctuations.


convolution of a real and complex signal How do we define convolution of: $$x(t)= \mathrm{sinc}(W t)$$ and $$y(t)= -j\mathrm{sinc}(2W t)$$

Answer: In the frequency domain, both of them are rectangular functions and multiplication of them is equal to the rectangular signal with smaller bandwidth i.e $x(t)$, which means convoluted signal is $\mathrm{sinc}(W t)$ with some scaling factor.

But I do not understand how this could be done, the two signals are orthogonal (because one is real and the other imaginary). So shouldn't convolution be 0? I would appreciate if someone could explain.

You are confusing several different notions. First of all, two signals $x(t)$ and $y(t)$ are orthogonal if

$$\int_{-\infty}^{\infty}x(t)y^*(t)dt=0\tag{1}$$

From (1) it is clear that if two real-valued functions $x(t)$ and $y(t)$ are not orthogonal, i.e. the integral in (1) is not equal to zero, then also $x(t)$ and $jy(t)$ are not orthogonal. So your argument that the two given signals must be orthogonal because one of them is real-valued and the other one is imaginary is wrong.

Furthermore, even if two signals are orthogonal to each other, their convolution is usually non-zero. Just take two rectangular functions that do not overlap. Obviously, they are orthogonal, but their convolution is a triangular function.

Consequently, one way to interpret the convolution of the two sinc functions is that a low pass signal with cut-off frequency $W$ is filtered by a low pass filter with cut-off frequency $2W$, and multiplied by the (irrelevant) factor $-j$. Since the low pass filter has no effect on the low pass signal, the only change to the signal is the constant factor $-j$.


frequencies in sound: multiple possibilities? First, I am by no mean a sound engineer (as you will guess later).

I was just wondering something while looking at the waveform of a .wav

for a given shape of waveform on a duration of 2 sec for example, how can we make sure that the frequencies fft gives are the only correct one?
what if very little parts of a sinusoid could be considered instead of a full continuous sinusoidal movement that just varies in amplitude?
or, what if a sinusoid had a lot of very fast varying amplitudes ?

that would lead to an infinity of solutions I guess..

The Discrete Fourier Transform, or DFT (the FFT is an algorithm that computes the DFT) of a length of finite duration (which any practical transform would need to be) is identical to the result of the transform of an infinitely long sequence formed by repeating the original sequence in time.

Knowing this provides the following insights related to your question:

First, anything that repeats (exact repetition) in time can only exist at discrete frequencies in the frequency domain. We see this with the Fourier Series Expansion specifically as shown in the graphic below.

The concept behind the Fourier Series Expansion is that any single valued continuous function can be represented as a sum of sinusoidal components and notably each component MUST have a frequency that is an integer of 1/T where T is the duration of the signal in time. Therefore IF those sinusoidal components were allowed to play out for all time (rather than being bound to the time interval [0,T]), the next cycle immediately after T would have to commence and proceed exactly as the waveform did at the start of the sequence (as each sinusoidal component would do the same). Thus it is often described that the Fourier Series Expansion decomposes any periodic function into a sum of sines and cosines (or equivalently and I believe mathematically simpler, complex exponentials).

The mathematical model of the time limited signal from [0,N-1] (for N samples of the DFT, which is equivalent to the analog bound of [0,T] in time), to also be a signal extending to infinity repeating with time provides for further intuition into the behavior and result of the DFT. Specifically we see that repeating a signal in time results in discrete uniform spaced impulses in the frequency domain, that can only exist at integer multiples of 1/T (including 0) where T is the length of the base waveform in time.

The second characteristic of the DFT is that it is done on a waveform that is sampled in time.  Without going into significant detail, sampling in time is associated with repetition in frequency (for those familiar with A/D and D/A conversion this will be readily apparent). So with the DFT we have both characteristics of repetition in time and sampling in time, which therefore means we will have "sampling" in frequency (discrete frequencies where only non-zero values can exist) and repetition in frequency.

The repetition is an implied construct that actually helps significantly to provide an intuitive understanding of many signal processing constructs- especially when considering both analog and digital domains. To be clear, the DFT involves only a fixed duration sequence both in the time and frequency domains, but mathematically these sequences can be repeated for infinite duration with the same result.

For example, to your question of what happens where there is a partial cycle? If we realize this equally represents a waveform repeating for all time, we see that in this case there would be an abrupt transition in the waveform.

Such a waveform cannot be created or represented with a single sinusoidal tone. Going back to the Fourier Series Expansion, we know that it can be represented by multiple tones as long as they are spaced at integer multiples of the repetition rate. So where we thought one tone exists, similarly in the DFT there will be multiple tones as required to create such an abrupt transition.  According to the DFT these frequencies really exist as what we are solving for in that process is the frequency components that are needed to create the time limited time domain waveform.

Below shows the case above (lower plot) compared to what we would get if there was a complete integer number of cycles over the time duration used (upper plot).  This is one explanation of "Spectral Leakage" with the DFT but given as example insight into the relationship between the time duration of the DFT chosen and the frequencies that would result.

If the waveform was changing with time (either in frequency, or amplitude beyond the sinusoidal component itself, such as an envelope) this would require many frequency components to represent it. This is no different than a modulation view of the time domain waveform: to transmit signals over the air we modulate the amplitude or frequency of a carrier frequency (that is better suited to go over the air) which results in several frequencies being present around that carrier. If the waveform is not actually repeating with time (which is likely the case), then instead of discrete tones we will get a continuous band of frequencies around the carrier.


Matlab filter design,through group delay I have the measurements of a filter’s Group delay and S-parameters.
The S-parameters are of the following form, presented in a touchstone file.

  !Agilent Technologies,N5242A,MY49421489,A.09.33.09
  !Agilent N5242A: A.09.33.09
  Date: Thursday, March 22, 2012 19:50:06
  !Correction: S11(Full 2 Port(1,2))
  !S21(Full 2 Port(1,2))
  !S12(Full 2 Port(1,2))
  !S22(Full 2 Port(1,2))
  !S2P File: Measurements: S11, S21, S12, S22:
   Hz S  dB   R 50
  1450000000 -0.44925556 132.79056 -43.664959 42.970737 -43.609634 45.291161 -0.41757283 131.60133

The filter’s Group delay data are like those demonstrated below:

  !CSV A.01.01
  !Date: Thursday  March 22
  !Source: Standard
  BEGIN CH1_DATA
  Freq(Hz)    S21 Delay(s)
  1.45E+09    -2.02E-08
  1.45E+09    -1.32E-08
  1.45E+09    -1.77E-08
  1.45E+09    -1.70E-08
  1.45E+09    -1.56E-08
  1.45E+09    -1.36E-08
  1.45E+09    -1.20E-08

What I want to do is simulate the impact that the group delay will have on a particular waveform using MATLAB. I tried to use the fdesign.arbgrpdelay in order to insert my group delay data and somehow observe how that would impact on a waveform but I am getting the following errors:

Error using fdesign.abstracttype/superdesign (line 96)
Design options must be specified as a structure or as parameter-value
pairs.

Error in fdesign.abstracttype/design (line 11)
    varargout{1} = superdesign(this, varargin{:});

Error in allpassfilterarbitrarygrpdly (line 413)
 Hgd = design(hgd,'iirlpnorm','Weights','MaxPoleRadius',0.95);

I also tried to use fdatool but I couldn't find a way of designing a filter by changing its group delay. The group delay was flat in all the available designs.

Does anyone knows how I can use the Group delay measurements in matlab to simulate a filter through its group delay?

The Solution finally was to match the input vectors in the group delay simulating Matlab functions.That means that the frequencies and the group delay Vectors given as input in Matlab should be SHORT and not as long as i had to give,and exactly on the same length (which makes sense).The most tricky part is to find which length was that would fit in and run the simulation.


Detecting sound inside a sound I was searching on the Internet hopelessly about some materials regarding detecting a sound in another sound. Say, I've got a recorded short sound (which may be anything, fragment of speech, fragment of music, noise, their combination or something else) and a longer sample, which contains this sound somewhere. The problem is, that the sound may be partially mixed with some ambient sounds. I know though, that despite mixing with other sounds, this will be exactly the sample I have recorded earlier.

Which algorithm may I use to perform such search?

I don't know if you are familiar with statistical signal processing and therefore will write my answer assuming that you are not. Everything I explain here is much better presented in any book about statistics. I would recommend Kay's book about detection theory. I first summarize your question by reformulating the 2 points you made, first in comprehensive English and then from a signal processing point of view.

You have access to two signals:

The signal you are looking for. Be it music, speech, or whatever. The only important assumption about this signal is that you perfectly know it. We will name it the signal of interest and denote it $s$. Let's say that it consists in $N$ samples denoted $s[n]$, with $0 \le n < N$.
A recording that contains the signal of interest as well as many other sources. The recording will thereafter be named the observation signal and denoted $x$. The other sources are all considered as noise in that they consist in anything but the signal of interest, and we will denote the noise $w$. The observation is longer than the signal of interest. Let's say it consists in $K$ samples denoted $x[k]$, with $K \ge N$ and $0 \le k < K$. Same goes for $w$.

You said you know what you are looking for ($s$), but you don't know  how badly it is "drown" ($w$) in what you are observing ($x$). From this point, the first thing I would try is a matched filter (MF) approach. I can't say for sure it will work for your problem because you gave us nothing about $w$. What I can say for sure is that this is the easiest approach that is very likely to work, both from a theoretical and a programming point of view. If this does not work, then you can come back here again. The fun part will start from there :)

But for now, let's focus on the matched filter. $s$ is shorter that $x$, so you should be able to take $M$ successive snapshots of $x$ of the same length as $s$, each denoted $x_m$. From $s$ and those $M$ snapshots, you can construct a third signal $y$ of length $M$, with each sample defined as follows :
$$y[m] = \sum_{n=0}^{N-1}{x_m[n]s[n]}$$

If $s$ is really in $x$, then you may see a clear maximum in $y$, at the index $m$ of the snapshot that contains $s$. If not, two possible conclusions:

$s$ is not in $x$. But that should not be the case, according to the presentation you gave of your problem.
Either the "level" of $s$ is way too low as compared to $w$ (we talk about a bad signal-to-noise ratio), or the statistical properties of $w$ are making this more complicated. The fun starts in the latter case.


What is the difference between complex FFT of baseband fft(I+jY) vs. FFT of absolute fft(abs|I+jY|)? If I demodulate an input signal at a frequency and receive its I and Y (in-phase and out of phase) values; I can also take the R value by finding the absolute of I and Y.
But what I'm confused is the application of FFT of a complex time signal (I and Y in this case) vs. FFT of a real value (R). I understand both result in a complex result, and with the FFT(R) has symmetry at 0 Hz between -BW and BW. And FFT(X+iY) does not have these symmetric properties anymore?
If I plot the absolute of FFT(X+iY), there would be no symmetry between -BW and BW, but what more information do I gain/lose from either methods?

Well, taking the absolute value is a nonlinear operation you do on your signal.
That means it introduces frequencies into your signal (harmonics!) that weren't there in the original, and you have no way of figuring out whether the things you observe are from your original signal or from your nonlinear operation.
So, you basically make your whole thing useless for spectral estimation. Don't do that!
It's easy to imagine how that goes wrong. Take the simplest to analyze signal we can think of: a simple real-valued cosine $\cos(2\pi t)$. We know the spectrum of that – a line spectrum, with components at $\pm 1$, each having power $\frac 12$.
If you take the absolute of that, that's (just for this specific signal!) equivalent with multiplying with a $\pm 1$ square wave that is scaled and shifted exactly so that it becomes negative whenever the cosine is negative. That way, the product is always the absolute value.
Multiplication of two signals means in spectrum a convolution of the spectra. The cosine spectrum we know, $\frac12(\delta(f-1)+\delta(f+1))$, and the spectrum of a square wave is an infinite sum of such deltas, at every odd multiple of the fundamental frequency $f$.
So, your "innocent" absolute just converted your pure and easy cosine to something with infinitely many spectral components. No good.


How do I interpret the result of a Fourier Transform? For example, I entered the following "equation" into Wolfram|Alpha:

FourierTransform[Piecewise[{{sin[t],t > 0 and t < 2*pi}}, 0], t, \[Omega]]

So as to take the Fourier Transform of a signal like this:

And turn it into some graphs that look like this:

My understanding from constantly — but casually! — hearing about Fourier Transforms (and DFT/FFT) in the context of audio processing and SDR I/Q data processing, is that I now have some information about a portion of the spectrum of the input signal. I'm trying to deepen my understanding, however, and refresh my memory of all I started learning in a Differential Equations class years ago.

Given my more recent practical background I'm left with the following questions understanding what the Fourier transform just did for me:

Why did I get a complex (real/imaginary) output from a real-only input? (How would I have passed some sort of "I/Q" equivalent in to the "equation" I provided?)
What units are omega, is it essentially Hertz if t is in seconds? What does negative Hertz mean?
If the input function were actually a voltage on a wire, would the real and imaginary parts be amplitude and phase (respectively) of the component frequencies?

From https://dsp.stackexchange.com/a/24758/18819, I gather that if I primarily care about power at a given frequency I would take the magnitude of the complex vector. From a Wikipedia article on negative frequency I'm not sure what its "wheel spinning backwards" analogy would mean in the context of e.g. electromagnetic or acoustic waves, but I'll take it to mean maybe it's just "some more of the signals" and combine the values across the y-axis.

So am I correct to interpret this as follows:

there is no power [correct term?] at 0Hz
at +0.5Hz and at +1Hz there is 1 unit (?) of "power" but
…at -0.5Hz there is also 1 unit of power
…while at -1Hz there is -1 unit of power
so at "real world" 0.5Hz there might be 2 units of power
but at "real world" 1Hz there might be 0 units of actual power?

Or what would be the best way to understand the first Fourier plot above, and derive relevant results from it?

The Fourier transform of a function $x(t)$ is defined by

$$X(\omega)=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt\tag{1}$$

Since $e^{-j\omega t}$ is complex-valued, the Fourier transform $X(\omega)$ is generally a complex-valued function, even if $x(t)$ is real-valued.

$\omega$ is the frequency in radians per second (rad/s), and it is related to the frequency in Hertz $f$ by $\omega=2\pi f$.

Since $X(\omega)$ is generally complex-valued, it can be represented by its real and imaginary parts

$$X(\omega)=X_R(\omega)+jX_I(\omega)\tag{2}$$

or by its magnitude and phase

$$X(\omega)=|X(\omega)|e^{j\phi(\omega)}\tag{3}$$

Clearly, the real and imaginary parts are not the same as magnitude and phase.

Note that for real-valued $x(t)$, the real part $X_R(\omega)$ and the imaginary part $X_I(\omega)$ of the Fourier transform $X(\omega)$ are even and odd functions, respectively. Consequently, the negative frequency components are redundant, and knowledge of the spectrum at positive frequencies is sufficient. If you are interested in power, you should rather plot $|X(\omega)|^2=X_R^2(\omega)+X_I^2(\omega)$ (which, by definition, is non-negative, as it should be the case for energy or power). The function $|X(\omega)|^2$ shows the distribution of the frequency components of $x(t)$.

As you've noted, the DC component of the signal is zero, because its average value over time is zero:

$$X(0)=\int_{-\infty}^{\infty}x(t)dt=0\tag{4}$$


Improving the performance of the Hough transform I have just started doing image processing, and I am currently working on line detection. I am using Hough transform on the image below.

I'm wondering how I can improve the efficiency of the Hough transform? Before applying the transform what can be done (pre-processes) to improve the performance?

Are you using openCV?
Anyway, if you are running an edge-detection, and then the hough transform on the edges, it might help to smooth/blur the edges with some bluring functions. OpenCV got some different ones for that.

In particular a boxFilter with a 2x2 sized kernel, anchor=Point(-1,-1), and normalize=false.
It will hopefully make all the edges "thicker", filling in gaps and inconsitencies in the process.


Gaussian weighted distance between pixel amplitudes: motivations and sources? In several image processing papers, and for instance in Random Walks for Image Segmentation (2006) by Leo Grady, changes in image intensities are turned into edge weights. In some fields, it is called Geographically Weighted Regression. This applies to graphs too, and has become popular for instance in non-local means denoising algorithms. A popular weighting is Gaussian. Quoting Grady (Equation 1, page 1772):

  In this work, we have preferred (for empirical reasons) the typical
  Gaussian weighting function given by:
  $$ w_{ij} =  \exp{\left(-\beta(g_i-g_j)^2\right)} $$
  where $g_i$ indicates the
  image intensity at pixel $i$

My questions are:

what are the earliest/most authoritative references for such a typical Gaussian weighting?
since Gaussian models are not considered the most appropriate for structured image modelling, what are the non-trivial "empirical" benefits of such a weighting, with respect to, eg, a near-Gaussian bisquare function?
$$w_{ij} = (1-\beta(g_i-g_j)^2)^2$$

Maybe not an exact answer, but I'll to give a direction.

What you are using is essentially an RBF-Kernel. First, it has a ready interpretation as a similarity measure (satisfies Mercer's conditions). In essence, other such kernels can be used. For instance, an inner product (linear kernel) would define an angular distance between pixels/features. For a deeper understanding, I would dive in the literature of kernels.

Having said that, the simplest form of such relations are called inverse distance weighting (Grady could have used this as well), where any identity like $1/d(g_i, g_j)^p$ can be used. These methods also cover a wide spectrum.


Correlation : How is correlation analogous to dot product? Till now I know correlation tells about similarity. I was watching a video lecture on image similarity in which I came to know that correlation is analogous to dot product. And hence correlation of two images is maximum when these images are similar as happens in dot product of two aligned(similar) vectors.

Dot product of two vectors $a=[a_1, a_2, a_3, .....]$ and $b=[b_1, b_2, b_3, .....]$ is given as $$a\cdot b=\sum_{i}a_ib_i$$ While cross correlation of two discrete real signals $x(n)$ and $y(n)$ is given as $$R_{xy}(n)=\sum_mx(m)y(n+m)$$
So now I want to know:

How correlation is analogous to dot product?
When is correlation maximum?
Why?

So using your formula for the cross-correlation:

$$R_{xy}(n)=\sum_m x(m)y(n+m)$$

we see that

$$\begin{align}
R_{xy}[n]&=\sum_m x[m]y[n+m]\\
&= \sum_m x[m] \tilde y_n[m]\\
&= x \cdot \tilde y_n\\
\text{for}\\
\tilde y_n[k] &= y[n+k]\text{ ,}\\
\end{align}$$
which is a dot product.

So for any "shift" in the cross-correlation, you have the dot product of the first operand with the $n$-shifted second operand.

Note, by the way, that your formulas are only correct for real-valued $x$, i.e. in general, the correct formula for dot product and cross correlation are:

$$\begin{align}
x\cdot y &= \sum\limits_{i=0}^{N-1} x^*[i]y[i]\\
R_{xy}[n] &= \sum\limits_{i=0}^{N-1} x^*[i]y[n+i]\\
\text{with}\\
x,y &\in \mathbb C^N\\
{z}^*&= \Re\{z\} - \Im\{z\}
\end{align}$$


Guitar Tuner or Guitar effects Algorithms? I need to explore the algorithm of the digital guitar tuner or the digital guitar effects algorithms.
Is there any Book or articles that you recommend of?

All I know is that relate to DSP.

thanks

Best book on the topic is probably http://www.amazon.com/DAFX-Digital-Udo-ouml-lzer/dp/0470665998.

There is some freeware stuff (very mixed quality) here: http://www.musicdsp.org/

Many companies spend a lot of effort in making good sounding effects and well working tuners to make money with it, so access to "good" non-commercial sources if fairly limited.


DTFT frequency range $$X(e^{j\omega}) = \sum_{n=-\infty}^\infty x[n] e^{-j\omega n} $$

The frequency term $\omega$ in DTFT is normalized as $\omega = \frac{\Omega}{f_\mathrm{s}}$

$\Omega= 2 \pi f$ is the angular frequency continuous-time Fourier Transform of $x(t)$.

A continuous-time signal that can be properly sampled will have non-zero energy components in its spectrum having frequency $|f| \lt \frac{f_\mathrm{s}}{2}$

Hence the range of $|\omega|$ should lie between $0$ through $\pi$. Given that the spectrum of a discrete time signal repeats every $2\pi$ radians, we should see this spectrum replicated from 0 to $\pi$, $2\pi$ to $3\pi$, $4\pi$ to $5\pi$ and so on.

I can extrapolate this to negative factors of $2\pi$ as well. So the spectra can additionally replicate from $-2\pi$ to $=3\pi$, $-4\pi$ to $-5\pi$ and so on.

I have however read that the range of $\omega$ is from $-\pi$ to $+\pi$. Where does this come from?

The DTFT is independent of Nyquist sampling. The DTFT is applied to a sequence of numbers regardless of any analog waveform that may be sampled. Economic sequences are often inherently discrete time sequences.  The DTFT applies to aliased sequences as well.

The summation is in general over $n$ from $-\infty$ to $\infty$.  The summation limits can be truncated to the nonzero range of $x[n]$.

The nominal convention as presented in Oppenheim and Schaefer, which they introduce prior to their chapter on Nyquist sampling, is that the DTFT is a continuous valued, periodic function defined uniquely on $-\pi$ to $\pi$ which has a length of $2\pi$.

For real valued sequences, the interval $-\pi$ to $0$ is conjugate mirrored of the spectrum over $0$ to $\pi$.

In general for complex valued sequences, the mirrored conjugate symmetry can’t be assumed.

$0$ to $\pi$ is only one half of a length $2\pi$ period.

edit

$-\pi$ to 0 has an image at $\pi$ to $2\pi$ and $3\pi$ to $4\pi$ and ......

0 to $\pi$ has an image from $2\pi$ to $3\pi$ and $4\pi$ to $5\pi$ and ....

similarly for the negative frequencies.


What is an "innovation filter"? I'm a math postgrad student working through a paper on eigenvalue decompositions of matrices of FIR filters (used for stuff like total decorrelation, convolutive mixing, MIMO). Towards the beginning, when describing their signal model, they use the phrase:

  If the PSD of the $l^{th}$ source is generated by a stable and causal innovation filter $\text{F}_{l}(z)$, ...

What does this mean? I understand a stable filter is one whose coefficients $a[n]$ approach 0 as $n \rightarrow \infty$, and a causal filter is one which does not use values from the future (i.e. no negative time lags). I have read that the innovation in terms of stochastic processes is defined as the difference between the actual next value of a time series, and the value given by an optimal prediction based on the information already available.

I cannot, however, find any description of an innovation filter. I have a mathematics background rather than a DSP/engineering background, and so sometimes some of the terminology trips me up, which is why I was hoping DSP stack exhange could kindly help!

An innovations filter of a WSS process $X(t)$ is a causal and stable minimum phase filter that can be used to generate $X(t)$ from a white noise input $N(t)$:

$$X(t)=\int_0^{\infty}h(\tau)N(t-\tau)d\tau\tag{1}$$

where $h(t)$ is the impulse response of the innovations filter.

The (causal and stable) inverse filter of $h(t)$ is called the whitening filter of $X(t)$. Its response to the input $X(t)$ is white noise $N(t)$.

A process that can be represented by $(1)$ is called a regular process. Note that its power spectrum is given by

$$S(\omega)=|H(\omega)|^2\tag{2}$$

where $H(\omega)$ is the Fourier transform of the impulse response $h(t)$ of the innovations filter (assuming that $E\{N(t)N(t+\tau\}=c\delta(\tau)$ with $c=1$). The requirement of causality of the innovations filter implies that the power spectrum $S(\omega)$ of $X(t)$ satisfies the Paley-Wiener condition

$$\int_{-\infty}^{\infty}\frac{|\log S(\omega)|}{1+\omega^2}d\omega\lt\infty\tag{3}$$

This implies that $S(\omega)$ cannot contain spectral lines, and it cannot be band-limited. Unlike a singular process (consisting of spectral lines), a regular process cannot be parameterized by a finite set of random variables and it is not predictable, i.e., it is not completely determined in terms of its past.

Reference:
Probability, Random Variables, and Stochastic Processes, A.Papoulis, Athanasios 1965. McGraw-Hill.


FMCW Radar Working Logic Question In FMCW radar rather than sending a signal at a constant frequency we are sending a chirp which goes e.g. from 6 GHz to 6.1 GHz in 5 ms.

What I cannot fit in my mind is this (numbers are just for example):

The RF signal is moving at the speed of light so while we are changing the chirp frequency linearly e.g. we start from 6 GHz and waited like 500 us, then increased it to 6.01 (up to 6.1 GHz) and we already received the 6 GHz signal from receiver without completing 5 ms chirp period. So up to this point we have nothing to do with other frequencies of chirp (6.02, 6.03, ..., 6.1). Basically how IF shift in signal occuring. Is it occuring for every frequency part of the chirp signal (6 GHz for 500 us, 6.01 for 500 us ... up to 6.1 GGz) or is it occuring as a result of whole chirp package?

Shortly should i just simultaneously change chirp from 6 to 6.1 and 6.1 to 6.0 back and forth and at the same time A/D convert the upcoming signal and save it somewhere until it makes 5 ms then plot the fft of the whole 5ms block.

I hope I am clear. Please show me what I am missing here :)

LATER UPDATE: Everyone is trying to make me understand the theoretical background where i already know but thank you guys for your time.

I am teling this with examples so people who have same problem in future can understand. While sending a chirp we start sending with a fix freq. as a start point like 6ghz. While it is traveling to the target and coming back a little bit of time passed and during this time our TX chirp signals freq. changed a little bit linearly like 6+X GHz so multiplexer difference is X Hz. Therefore this change depend on the distance of the target since let say target is so far far away we will maybe receive the first signal which was at 6GHz while our chirp is at its halfway of 6 to 6.1 GHz cycle. So the difference will be like 6.05-6=0.5GHz.

Another way of saying, for stationary object nothing happening on the transmitted signal(doppler effect says reflected data is radited with same structure if object is stationary) but until the transmitted signal come we changed the transmitted signal so the upcoming signal will be multiplexed with different freq. signal and based on this change we are telling this is object at this much distance. Thank you

Grace Hopper used to use pieces of optical fibre, to show how "long" is a nanosecond at the speed of light. It is approximately 30cm and a nanosecond is a sort of a time interval that can be easily measured, accurately, without expensive lab equipment (e.g. with a simple counter, using a crystal as its time basis).

So, the "speed of light" is not something that is inaccessible or somewhat of a big problem for electronics. Especially at Radio Frequencies (RF), the speed of light is even a "necessity" if you are trying to observe phase shifts that imply extremely small distances in space.

The point about the continuous wave radar is that it can use its local oscillator for the "demodulation" of its transmitting wave. And using different patterns of modulation, it can return different types of information (for example range and / or velocity of a target).

So, at its heart, the Continuous Wave radar is based on measuring phase shift via a phase detector. And again, this operation, at its heart, is extremely simple: Multiply two sinusoids at frequencies $f_1,f_2$, you get a signal with two components, one at $f_1+f_2$ and one at $f_1-f_2$. At the special case where $f_1=f_2$, you get a component at $2 \cdot f_1$ and a component at 0 Hz, or a DC component, a constant value added to the signal. If, in addition, you vary the phase between $f_1,f_2$ (while still keeping $f_1=f_2$), then the DC component is proportional to the phase difference.

This last point is also the basic principle by which Phase Locked Loops produce the error signal that they use to tune their local oscillators to the frequency and phase of an incoming signal.

So, having said (and hopefully understood) all this, let's look at what you ask:

  Basically how IF shift in signal occuring.

The shift is occuring because of phase shift. Another way of viewing phase ($\phi$) is as the first derivative of frequency ($f$). Therefore, if you integrate phase, you get frequency. In other words, many many many little phase changes, accumulated in the right way ammount to whole changes in frequency. For an explanation of why does the shift occur in the physical world and within different propagation media, please see the Doppler Effect.

  Is it occuring for every frequency part of the chirp signal (6 GHz for 500 us, 6.01 for 500 us ... up to 6.1 GGz)...

Yes. The same phenomena that apply at some $f$, also apply for some $f + \Delta f$, where $ \Delta f$ is some frequency offset. So, whether looking at a single sinusoid or a set of sinusoids traveling together, there will be phase changes observed, proportional to each frequency's wavelength. Which brings us to the second part of this question:

  ...or is it occuring as a result of whole chirp package?

Again, yes, in the case of the Continuous Wave radar, as it was mentioned before, you use the local oscillator to demomdulate the incoming (reflected) wave. Therefore, your "window of integration" (the window within which you assess phase change) would have to be proportional to the "length of the pulse".

Very briefly: If you get very early reflections, at timings that would cause a negligible phase change between the transmitted component and its reflection then that would register a phase change of zero and consequently also a default distance. In other words, the frequency and pulse timing of the CW radar set limits to the distances it can observe. Any reflection by targets at LESS than that critical distance will not be registered and any reflection by targets at GREATER than that critical distance will also not be registered (or contribute to noise if they happen to be received within a subsequent "window"). It is reflections within specific limits that cause detectable phase changes which can then be used to infer either range (distance to target) and / or velocity of target with respect to that of the source (if they are both moving).

For more information, please see differences between "types" of CW, purely by changing the "pulse" pattern, via this link.

Hope this helps.


The phase of signal when it passes a frequency divider? When we used the microwave frequency divider, such like: RF Bay FBS-N-14 or analog dividers in wiki.

The input signal is $\cos(\omega t+\phi)$, and the output signal is $\cos(\frac{1}{2}\omega t+\theta)$ when it passes the divide-by-2 prescaler.

what is the phase $\theta$ equal to? $\frac{\phi}{2}$ or $\frac{\phi}{2}+\pi$?

That can't generally be said. It depends on the actual design of the device.

In fact, it's not going to be constant for all frequencies, because that's physically pretty much impossible to achieve.

You'll find multiple descriptions in literature for three- and four-port devices like these, and you'll have to learn how to infer the relative phases for different frequencies from the measurements that the manufacturer of a divider will provide to you.


What's the meaning of ergodicity? I just read the topic about Ergodicity but I have ambiguity about its meaning (by intuition). What does mean: (for mean) Statistical average = Time average. Could you please explain it in detail. Thanks.

If you sample a random process for a specific t, you will get one realization of a random variable. For another t, you get another realization of that random variable. This random variable has its statistics which is almost impossible to learn in real world because not all sample paths are observable. See the brown rectangle in the figure below.

That becomes possible in ergodic processes. An ergodic process is the one where the time average of a sample path is the same as the statistical mean. See the purple line in the figure below, which was the observed sample path such as a particular realization of noise in a communications receiver.


How to calculate Noise and distortion of an image? I need to make an application in iphone which needs to calculate noise, geometric deformation other distortions in an image. How to do this? I have done some image processing stuff with opencv + iphone. But I dont know how to calculate these parameters.

1) How to calculate noise in an image?

2) What is geometric deformation and how to calculate geometric deformation of an image?

3) Is geometric deformation and distortion are same parameters in terms of image filter? or any other distortions available to calculate an image is good quality or not?

You should do it on a known set of calibration targets. For example, this type of calibration target is being used to measure noise. You should find the patches locations, and measure the noise by some metric. If you don't know any metric, you can start by browsing DXOlabs site, or Imatest site.

I guess that you mean radial distortion. Radial distortion can be found by a process called camera calibration, by using a known target such as checkerboard. For more info, see Camera Calibration Toolbox. Here is an image from another software that does it:

I didn't understand your question at all. Please reprase it and add more information.


What is the role of Freqency in Digital Images Possible Duplicate:
  What does frequency domain denote in case of images?

If we look deeply into Digital Images they are purely generated by the pixels distributed in rows and coloumn and the pixel consists of number of certain bits which define the complete Image.

Then what is the role of Frequency in generating Digital Images ? like many people says that in Digital Image processing and Image Compression, Wavelets is much more efficient than Discrete Cosine Transform as Wavelets gives the information of both time and Frequency.

I want to know what role does frequency and time plays in producing a Digital image, does it has to do something with pixels?

You can take images as 2D discrete signals. The "time" in 1D signals is actually two spatial dimensions in images (2D signals).

You can measure "frequency" as well - imagine a line of white pixels with regular spacing. The spacing represents period $p$, and frequency is given by $1/p$. Hence the maximum frequency the discrete signal can contain is limited by the sampling rate.

When creating a digital image, sampling can cause new frequencies to appear due to aliasing. See Nyquist frequency.

Fourier and wavelet transforms have discrete variants to work with discrete signals - very simply said: integrals are replaced by sums.

Wavelets are more efficient in image compression, because they can deal with localized redundancies. Simply look of DFT of a grey square on black background and then on DWT of that square.

DWT will exploit transitions of the square shape in high fequency layer, and other layers will be almost clear. The sharp 1px transition would be best described with Haar wavelet without any remainder.

On the other hand, DFT or DCT will result in many cosine functions which summed together form the edge without visible ringing artifacts (Gibbs effect).

Wavelet have its artifacts too, but these are better localized, hence does not corrupt rest of the image (or rest of the macroblock in case of JPEG/DCT compression).


recognize unwanted sounds, noises, and deleting them algorithm I have a function that I need to use denoise algorithm on.

What I thought is to divide it to smaller parts and then use Fourier Transform on each one of the parts, but I am not sure what to do next.

Can you please help me and give me few algorithms to recognize noises, if possible algorithm that will fit the structure I described above?

Provided that the noise is uniform across your sound and that you have regions where only noise is present, you can make use of a technique called spectral subtraction. The whole process consist of 4 steps.

  Find a noise-only region and apply an FFT on it in order to obtain the noise profile (noise spectrum). The longer the noise region, the better the noise profile (resolution) you'll get (use longer FFT lengths)
  Apply FFT on your entire sound. This way, you'll get the spectrum for both the signal (people talking) as well as the noise.
  Subtract the noise profile from the overall sound spectrum (2-1). Ideally, you'll end up with the signal-only spectrum (of course, some
  noise will still remain).
  Perform an inverse FFT on the signal spectrum obtained in step #3 to get a time-domain signal (waveform) .

If the filtered signal still has some unwanted artifacts, you may try to interpolate the signal spectrum prior to performing the IFFT (step 4)

Note that you don't have to break up your signal into smaller chunks (this is not STFT)


Power of summed signals Pretty much a newbie here, but would like to understand a pretty basic concept.

We all know that a quick and dirty way of calculating signal power of a part of a discrete sampled system (such as audio signals) is the age-old power formula for RMS power, where we simply break the signal into chunks of N samples and sum the squares of the amplitudes, then divide by the number of samples

$$
P_\mathrm{x} = \dfrac{1}{N}\sum_{n=1}^{N}|x[n]|^2
$$

So let's say I have 10ms samples of two completely uncorrelated signals at 48khz (480 samples), and I have calculated the RMS power of these samples using the above formula and have come up with P1 and P2, the power of my small segment of each signal.

Then if I sum the two signals (simple software mixer function), what is the resultant power of the two summed signals?  Is it simply P1 + P2?

If so, do I simply keep going?  So if I had four signals with P1, P2, P3 and P4, if I sum all four of them together, is the resultant RMS power P1+P2+P3+P4?

The formula you gave is an estimator for the "statistical power", i.e. the expectation of the square of a continuous  underlying random process:
$$ P_X = E(X^2(t)) $$
Computing the latter theoretically would require to know the probabilistic behavior of $x(t)$, which you don't have access to in real life, especially since:

you get sampled data (not continuous),
you generally have a limited number of realizations.

This is even more complicated when you add variables, because you should know their paired probability behavior. A classical result is that, for variables $X$ and $Y$,
$$\mathrm{Var}(aX+bY) = a^2\mathrm{Var}(X)+ b^2\mathrm{Var}(Y)+ 2.a.b.\mathrm{Cov}(X,Y)$$
This quantity can take different values, depending on the last term: $\mathrm{Cov}(X,Y)$. For uncorrelated variables, the covariance vanishes, and the variables add as:
$$\mathrm{Var}(aX+bY) = a^2\mathrm{Var}(X)+ b^2\mathrm{Var}(Y)$$
So, for uncorrelated signals, powers add. And this generalizes to several signals.
Note that when  you sum over a time frame to estimate such quantities, you somehow assumes that you can estimate quantities by averaging in time, with assume that the signals are correctly sampled for a stationary and ergodic process.
For $K$ signals of  identical power $P$, the combined power can take value:

minimum: $0$ when signals sum to zero,
$KP$ when uncorrelated,
maximum: $K^2 P$ when signals are identical.

Additional pointers:

Addition of signal variance and noise variance
Variance of Uncorrelated Variables
Signal averaging
Superposing Signals


Damped Harmonic Oscillation as an LTI The goal is to create an LTI filter which is exactly, or approximates, damping of harmonic modes.

The equation of course is:

$$\frac{d^2 x}{dt^2} + 2 \xi \omega \frac{dx}{dt}+\omega^2x=0$$

This can be done with very costly convolution, I suppose.

I was thinking that instead of treating the velocity term as a part of a time dependent potential (no longer time invariant), I can treat this as a two dimentional state and therefore eliminate second derivatives and have it a $z^{-1}$ feedback filter. It is working weird.

A continous alternative of something close is:

$$e^{-\xi t} \cos \omega t$$

How is this approximated efficiently?

Links are enough.

This is a classic dynamics problem. I think you do not have to go through the convolution to solve these type of problem. Instead, you can draw the frequency response function (FRF) of the system and multiply with the FRF of the forcing function to obtain the resulting FRF. Then the resulting FRF can be transformed into the time domain using the inverse Fourier transform. You can also solve it by doing the convolution. Then you have to use the Duhamel's integral to solve this type of problem.

The link is attached.

http://people.duke.edu/~hpgavin/cee541/sdof-dyn.pdf


Multitone signals with low crest factor I am trying to generate a multitone signals with low crest factor and as suggested in the Boyd's paper I am using the quadratic phase distribution:

$Phase(k) = \pi(k - 1)(k - 1)/N$

where $N$ is the number of tones in the signal and $k$ their frequency indexes.

This formula works well, giving the expected crest factor ($CF <= 2$), only if the frequency indexes are strictly sequential (e.g. $k = 2,3,4,...$) but if the indexes are not sequential (e.g. $2,5,7,11,...$) or, even worst, logarithmically distributed (e.g. $1,2,4,8,...$) the resulting $CF$ is very often much higher than $2$.

I played quite a lot with the above formula, trying differ combination of $k$ and $ N$, but without any success.

What is strange, to me, is that the lowest $CF$ is always obtained by using sequential $k$ values ($k = 1,2,3,4,...$), regardless of the actual position of the frequency index.

Just to complete the story, using the Shapiro-Rudin sequence I get more or less the same behavior.

Now my question is: do you know a way to amend the above formula to get low $CF$ with any arbitrary frequency index distribution or another phase generator algorithm suitable for this scope?

Ciao and many thanks for the attention.

Franco

okay, there is a slightly newer paper by Ivo Mateljan titled Signal selection for the room acoustics measurement.  i ended up reviewing the paper for the 1999 IEEE Mohonk conference, and when it turned up that Ivo was not going to make it to the conference to present the paper, they ask a reviewer to sorta stumble through the issues of the paper.  i ended up talking about what's wrong with MLS more than about Ivo's technique.

if you cannot get a copy of Ivo's paper (it costs money from IEEE), email me and i will send a pdf to you.  in that paper he has a section titled The Low-Crest Multisine which i paste in the text below.  perhaps you can modify it to fit with your selected frequencies.  i dunno.

i still haven't figured out SE's formatting to make the code look like code.  someone else can edit this answer to do that.

Various algorithms were proposed for lowering the multisine crest factor. The following algorithm is slight modification of algorithms proposed in [9] and [10]. It has fast convergence and resulting multisine has good phase randomisation property.

Algorithm: Low-crest multisine (LCMS) generation.

Generate uniform random sequence u of length N.
Transform the sequence to random sign sequence s:

    for i = 1 : N

        if (u(i) > 0.5)

            s(i) = 1

        else

            s(i) = -1

    end

Generate the multisine using the phase values of DFT(s).
If the number of iterations is less then 10, clip the time series to 85% of the maximum value, else clip to 90% of the maximum value.
Apply DFT to the clipped time series.
Generate the new multisine using the phase values obtained in the step 5.
Calculate the crest factor. If it is higher than objective value go to the step 4.

Although there is no proof for the convergence of this
algorithm, it has never failed in a more than thousand runs with
different starting seed factor. The number of iterations grows
exponentially when decreasing crest factor. For a multisine
sequence of length 32768 it needs only 4-5 iterations to halve the
crest factor from 4 to 2, but to get the crest factor < 1.5 it needs
several hundreds or more than thousand iterations.

Phase randomisation properties of this signal can be
estimated by testing how it distributes nonlinear distortion in the
estimated PIR. The following numerical experiment is done.
The input sequence of length 4096 drives 9th order bandpass
Butheworth filter (0.1-0.3 fs). The output signal is distorted with
2nd and 3rd order nonlinearities, which are equivalent to 1% of 2nd
harmonic distortion and 0.7% of the 3rd harmonic distortion.

Impulse response is estimated using crosscorrelation and scaled
to maximum value 1. Tables I and II show the rms and peak level
of the impulse response distortion (in percentage of the maximum
value). First 50 values of the impulse response are excluded from distortion calculation to better estimate peak/rms tail distortion. Also, the crest factor of the generated discrete sequence and a crest factor at the filter output are shown.

These results show that low-crest multisine generates the same level and a distribution of distortions as a true random phase multisine, regardless the crest factor value. MLS sequence generates slightly lower distortion rms level (1-1.5dB), but peak/rms ratio is much higher, especially with 2nd order distortion (10dB higher than with multisine excitation). This high peak/rms ratio disqualifies MLS as a reliable excitation signal for echo detection.

For all type of an excitation signal the crest factor at the filter output is almost the same. This is a general conclusion for narrow band systems. The benefit of using signals with lower crest factor is significant only for wideband systems, and when it is necessary to use all of the available measurement dynamic range.

RPMS is the only signal, which do not exhibit a significant crest factor transformation. This fact, constant amplitude and true random phase property fully qualify it as a periodic white noise signal.

[9] Schoukens, J., Pinelton, R., Ven der Ouderaa, E., and Renneboog, J., "Survey of Excitation Signals for FFT Based Signal Analysers", IEEE Trans. Instrumentation and Measurement, vol. 37, September 1988.
[10] Ven der Ouderaa, E., Schoukens, J., and Renneboog, J., "Peak Factor Minimisation of Input and Output Signals of Linear Systems", IEEE Trans. Instrumentation and Measurement, vol. 37, June 1988.


Understanding the definition of signal's bandwidth In class, we have defined bandwidth of a signal as the difference between the maximal and minimal frequency. Also, our tutor explained that the bandwidth indicates how fast the signal is changing.

We have been given an example of a single sinusoid which has zero bandwidth (since it's both the maximum and minimum frequency)

What isn't clear to me is why a rectangular shaped signal has infinite bandwidth.

How can it be explained rigorously and intuitively?

"tutor explained that the bandwidth indicates how fast the signal is changing"

This is not correct even in a hand-waving fashion.  As you were told, a pure sinusoid $x(t) = \sin(2\pi f_0t + \theta)$ has zero bandwidth (according to the definition of bandwidth that your tutor is using) but that does not mean that the sinusoidal is not changing at all! It is changing and its maximum slope (a.k.a. "how fast the signal changing") is
$$\max \left| \frac{\mathrm dx}{\mathrm dt}\right| = \max \left|2\pi f_0 \cos(2\pi f_0t+\theta)\right| = 2\pi f_0.$$ More generally, a result due to Bernstein says that if the maximum frequency in a generic $x(t)$ bounded within $[-1,1]$ is $f_0$, that is, $X(f) = 0$ for $|f| > f_0$, then
$$\max \left| \frac{\mathrm dx}{\mathrm dt}\right| \leq 2\pi f_0.$$
In short, the minimum frequency (in a bandpass signal) is irrelevant; it is the maximum frequency that matters, and the "how fast the signal is changing" depends on the maximum frequency and not on the bandwidth.

Turning to the rectangular pulse that puzzles you, note that the signal doesn't have a derivative (as far as mathematicians, including Bernstein, are concerned) when the signal makes the abrupt transition from high to low or vice versa but here on dsp.SE we have no qualms about saying that the derivative is $\pm \infty$ at these points.  As MattL says in a comment, if you know the Fourier transform of a rectangular pulse (hint: $\text{rect}(t) \leftrightarrow \text{sinc}(f)$), then you know that the rectangular pulse has no maximum frequency (there is no finite $f_0$ such that $\text{sinc}(f) = 0$ for all $|f|>f_0$) and that the infinitely fast changing of the rectangular pulse is well-supported by the infinite bandwidth of the rectangular pulse.


Coupled system instability Suppose there are two linear devices (containing some multi-physics but the details don't matter here), and these devices are given to us as two black boxes that can be studied experimentally. Black box #1 takes signal X as input and produces signal Y as output. Black box #2 takes signal Y as input and produces signal X as output. It is conceptually similar to a microphone converting an acoustic signal to an electric signal, and a loudspeaker converting an electric signal to an acoustic signal. If these two black boxes get connected to each other then an instability can probably arise. The question is what measurements on the two black boxes would produce all necessary and sufficient information to assess whether the coupled system would go unstable, and at what frequency. Would the right approach be to study the response of each black box to a sinusoidal input? Or the response to a delta-pulse input? What is the criterion for instability for the coupled system, and what measurements and calculations would be needed for assessing it?

Assuming both black boxes are LTI systems, that's reasonably easy to do.
Let's call the transfer function of the first box $H_A(\omega)$ and the second one $H_B(\omega)$. The open look transfer function is simply the product of the, i.e. $H_O(\omega) = H_A(\omega) \cdot H_B(\omega)$
A sufficient stability condition is $|H_O(\omega)| < 1$. For an actual acoustic feedback loop that's really all there is to it since there is always plenty of group delay in the system and the phase is spinning so fast that any type of phase criteria is pointless.
So measure the transfer functions (fully calibrated !), multiply them and check whether the max gain is smaller than unity.


How to implement an RF Mixer for chirp signal in python? According to this link, all I need to do to mix two signals is to multiply them.
I have two chirp signals with same properties, one leading the other. If I multiply them to get the IF output, I don't get a constant frequency which is what I expect, I get double the frequency:

from scipy.signal import chirp, spectrogram
import matplotlib.pyplot as plt
import numpy as np
import math

def plot_spectrogram(title, w, fs):
  ff, tt, Sxx = spectrogram(w, fs=fs)
  plt.pcolormesh(tt, ff, Sxx, cmap='gray_r')
  plt.title(title)
  plt.xlabel('Second(s)')
  plt.ylabel('Hz')
  plt.grid()

fig, ax = plt.subplots()

T = 0.4
fs = 7200
t = np.arange(0, int(T*fs)) / fs
f0, f1 = 5, 120
sc = chirp(t, f0=f0, f1=f1, t1=T)
rc = chirp(t, f0=f0, f1=f1, t1=T, phi=90.)

plt.figure(1)
ax.plot(t, sc, 'blue', lw=1.25, label='sent')
ax.plot(t, rc, 'green', lw=1.25, label='recv')

plt.figure(2)
plot_spectrogram(f'Chirp mul, f(0)={f0}, f({T})={f1}', sc * rc, fs)

plt.show()

Wouldn't mixing them give you a sum and difference frequency and you have to lowpass filter to get the difference frequency only?  Also your chirps are the same frequency, just out of phase?  So the difference would be zero and the sum would be double.
So:

Actually create a chirp that is a different frequency from the original, not just out of phase.
After mixing, lowpass filter to keep only the difference frequency


How to analyse periodic data obtained from a vibration experiment? I have data from a vibration experiment on a structure. This data was obtained by attaching an accelerometer to a FFT analyser, which records the data at a high sampling rate (typically in kHz range) and can also perform real-time FFT of the recorded data.

The image above shows the recorded signal below and the FFT above. The FFT shows all the dominant frequencies in the original signals.

What more patterns/ important data points can I discern from this data, that might have utility ?

For example, one being, the RMS value of the signal. What else ?

Also, the recorded data here as a mat file, which can be loaded in matlab, the data in a variable called Track1.

The physics of the problem: The data is collected from a shake table, a rigid platform, that oscillates back and forth, so the sharp edges in the original signal periodically, representing the point of return.

Remember that you are observing a system (the structure) driven by an input (shake table).  Assuming linearity, the FFT result was caused by the shake table and the structure.  You should be able to pick out the various components to identify mechanical resonances in the structure.  For instance, you should see a strong component at the frequency of the shake table.  Exactly how to do this in more detail is out of my domain, as I am an electrical, not mechanical, engineer.


Why complex signal has no imaginary spectrum I am learning about complex sampling.

I am confused why $~e^{ j 2\pi f~ n}~$ has only a real spectrum. I would have thought the $j ~\sin(2 \pi f n)$ would produce a single spike in imaginary spectrum just like there is a single spike in real axis from $\cos(2 \pi f n)$.

I understand that the spectrum is one sided because the negative complex exponentials cancel out, but why is there not a one sided real and imaginary spectrum?

Many thanks

The imaginary part of the spectrum corresponds to the odd part of the time-domain sequence. Since the given signal is even, i.e.,

$$x[n]=e^{jn\omega_0}=x^*[-n]=e^{-j(-n)\omega_0}=e^{jn\omega_0}$$

the imaginary part of the spectrum is zero, i.e., the spectrum is purely real-valued.

Note that for complex-valued signals, even in this context means $x[n]=x^*[-n]$, i.e., its real-part is even and its imaginary part is odd.

In sum, if a sequence $x[n]$ satisfies

$$x[n]=x^*[-n]$$

(i.e., the sequence is even and, consequently, its odd part is zero), then its DTFT is purely real-valued.


How to extract the maxima and minima of a noisy period input signal I have a (rough) known period sine wave coming in, held in an array, with superimposed noise. I need to find and count the maxima and minima (in C). Any suggestions or algorithms? The signal has already been received and buffered, so not real time.

If you have the entire signal, you can then run a Lowess smoothing. This algorithm is statistically robust and provides locally acceptable smoothing. Outliers will be tolerated as well. After smoothing your signal in that way, you might apply a derivative based (or even derivative free) peak detection in order to obtain the signal peaks. These will be your local minima/maxima.

You could use MATLAB or R and port smoothing algorithms to C, they are pretty straightforward and open source:

http://www.mathworks.com/help/curvefit/smooth.html

And here are your options in C:

1) https://code.google.com/p/variationtoolkit/ for Lowess/Loess.
2) For peak detection, there are many tools available in C, one of them being:

https://github.com/xuphys/peakdetect


What's the fastest way to get the phase angle and amplitude of a COSINUS 50hz signal? I'm not a DSP specialist.I'm working in transmission and distribution.
I need to retrieve the amplitude and phase angle of a discrete 50Hz COSINUS signal as fast as possible. Right now I'm using a Sliding DFT an I can manage to get this information in 10ms (half period) but I would like to know if it's possible to have it faster. I can adjust the sample time up to 1micro second.

I need to have this information (amplitude and phase angle) as fast as possible because it will determine the time constant of my corrector.

Thank you in advance for your precious answer

Depending on how the signal will look like you can get away with only two or three samples:

$$
y(t) = A + B \cos(\omega t + \theta)
$$

If you know what the main value, $A_0$, would be, then only two samples would be required. Namely the cosine, shifted by a phase angle $\theta$, can be rewritten as a linear combination of a unshifted sine and cosine:

$$
y(t) = A_0 + A_1 \cos(\omega t) + B_1 \cos(\omega t) = A_0 + \sqrt{A_1^2 + B_1^2} \cos\left(\omega t - \tan^{-1}{\frac{A_1}{B_1}}\right)
$$

The gains can be solved for by solving the following linear system of equations:

$$
\begin{bmatrix}
1 & \sin \theta_0 & \cos \theta_0 \\
1 & \sin \theta_1 & \cos \theta_1 \\
1 & \sin \theta_2 & \cos \theta_2
\end{bmatrix}
\begin{bmatrix}
A_0 \\
A_1 \\
B_1
\end{bmatrix}
=
\begin{bmatrix}
y_0 \\
y_1 \\
y_2
\end{bmatrix}
$$

where $\theta_i=2\pi f(t_0+i\Delta t)$, with $f=50\ Hz$ and the $t_0+i\Delta t$ are the consecutive moments in time used for for the three measurements and $y_i$ is the measured signal at those moments in time. The values found for $A_0$, $A_1$ and $B_1$ can then be used to find the constants of the initial equation, namely:

$$
A = A_0
$$

$$
B = \sqrt{A_1^2 + B_1^2}
$$

$$
\theta = -\tan^{-1}{\frac{A_1}{B_1}}
$$

This should work for any different points in time, which are not a multiple of the period of the signal spaced apart, however this does assume that there is no noise present in the signal. If there is noise then you might want to consider using more points and performing a least squares method. Also a shorter sample time/higher sample rate would make these calculations more sensitive to the noise, but you will be able to take more samples in the same amount of time.


Find the ber for the cdma protocol Do you know a reference (book, article) which gives (with a proof) the theoretical BER in the CDMA protocol for a differential modulation?

CDMA is not a protocol. It's a multi-access scheme.
A BER can't be generally given for a multi-access scheme. You'll need to define your access characteristics, SNR scenario, look at the specific multi-user detector for your specific modulation.

Hence, BlackMath is right: you'll need to go down the long route: Learn about CDMA's basics; do a lot of estimation theory, and apply it to the problem that you're looking at.

For many types of CDMA, detection probabilities are relatively easy to derive; things get a bit complicated because you'll need to define how the different users access the channel, in terms of occupation probability, and synchronity, as well as probability distribution of receiver powers.

Then you'd need to model how your multi-user detector deals with that. As you can see, this requires stochastic modelling of your scenario.


One integral inverse CWT MATLAB's icwt docs state inversion to be done by a single integral:
$$
f(t) = 2 \Re e\left\{
   \frac{1}{C_{\psi, \delta}} \int_0^\infty \left< f(t), \psi(t) \right> \frac{da}{a} \tag{1}
   \right\}
$$
Where does this expression come from, and how is integration over $b$ (translation) eliminated (as in Wiki's formula)? Does it hold if computing discretely over log-distributed $a$ (e.g. $(2^{1/32})^{[1, 2, ...]}$)?

The formula is premised on the wavelet being analytic, or being nonzero only over non-negative frequencies: ${\hat\psi} (\omega < 0) = 0$. (Note all wavelets also have ${\hat \psi (0)}=0$ per the admissibility criterion). From Daubechies et al,
$$
\begin{align}
\int_0^\infty W_f(a, b) a^{-1} da & = \frac{1}{2\pi} \int_{-\infty}^{\infty} \int_{0}^{\infty}{\hat f}(\xi) \overline{ {\hat \psi}(a\xi)} a^{-1} da d\xi  \tag{1a} \\
&= \frac{1}{2\pi} \int_0^\infty \int_0^\infty {\hat f}(\xi) \overline{ {\hat \psi} (a\xi)} e^{jb\xi} a^{-1} da d\xi  \tag{1b} \\
&= \int_0^\infty \overline{ {\hat \psi}(z)} \frac{dz}{z} \cdot \frac{1}{2\pi} \int_0^\infty {\hat f}(\xi) e^{jb\xi} d\xi  \tag{1c} \\
\end{align}
$$
where a rather confusing change of variables was replaced with a clearer one, $z=a\xi,\ dz=\xi da$. Note that this c.o.v. is only valid if integrals share limits, which is enabled by $\psi$'s analyticity; otherwise, $z$ would integrate $-\infty $ to $\infty$, but $z$ includes $a$, which only spans $0$ to $\infty$. An identity was also utilized in splitting integrals.
Assuming $f$ is real, we have ${\hat f(\xi)} = \overline{\hat f (-\xi)}$, thus from $0$ to $\infty$,
$$
\frac{1}{4\pi} \Re e \left\{ \int_0^\infty {\hat f(\xi)} e^{jb\xi} d\xi \right\} = f(b),
$$
and setting $C_\psi = \int_0^\infty \overline{ {\hat \psi (z)}} z^{-1} dz$ thusly derives
$$
\begin{align}
f(b) &= \frac{1}{2} \Re e \left\{ \frac{1}{C_\psi} \int_0^\infty W_f(a, b) a^{-1} da \right\}  \tag{2a} \\
     &= \frac{1}{2} \Re e \left\{ \frac{1}{C_\psi} \int_0^\infty f \star \overline{\psi_a(b)} \frac{da}{a} \right\}  \tag{2b} \\
     &= \boxed{ \frac{1}{2} \Re e \left\{ \frac{1}{C_\psi} \int_0^\infty \left< f(t), \psi(t) \right> \frac{da}{a} \right\} }  \tag{2c} \\
\end{align}
$$
Since $W_f$ sweeps $b$ over all $t$, above equivalently holds for $f(t)$.

The paper you linked uses $a^{-3/2}$ -- Indeed, as it assumes an L2 norm on the wavelet, whereas MATLAB uses L1; further info can be found here.
What about the $\delta$? -- The $\delta$ is from MATLAB writing the expression as a special case of the more general inverse, as it shows above near it, with $\left<f, g\right>$, which is derived here.
and the $2 \Re e$? -- Good question - perhaps MATLAB's $C_{\psi, \delta}$ is a bit different, but per above definition we clearly divide by two, not multiply. If in doubt, try both.

Interpretation: it helps to observe that $b$ in all derivations can be replaced with a constant - so the relation holds pointwise; if we integrate along $a$, i.e. a column of the 2D time-frequency representation, we recover the function at that point, $f(b)$.
Notice in the double-integral iCWT we use a dual wavelet to "undo" the inner product with the CWT wavelet, analogous to $e^{+j\omega t}$ in $\mathcal{F}^{-1}$; here, we don't use any inverting kernel - the transform undoes itself when integrated over all scales. This implies some symmetry property on the wavelet $\psi$ in encoding information across scales - though I'm unsure how this exactly works out; the most important is analyticity which permits this in the first place, so one may search there for an answer.

Discretized case with exponentially distributed $a$
From same page of linked paper, the discretized linear case reads
$$
f(b) \approx \frac{1}{2} \Re e \left\{ C_\psi^{-1} \sum_k W_f (a_k, b) \frac{(\Delta a)_k}{a_k} \right\},  \tag{3}
$$
with $(\Delta a)_k$ taking role of $da$. For the linear case, $(\Delta a)_k = \Delta a$, but not for the log case, as the increment size increases for greater $k$. That's the only distinction between log and linear.
In practice, however, it's better to directly use the discretized log-scaled transform - which derives as follows:
$$
\begin{align}
& z = e^a, \ dz = a da \\
& \frac{da}{a^2} = \frac{d \ln (z)}{z} \tag{4}
\end{align}
$$
Suppose $z(a) = 2^{a/32}$. Note that the log-plot of an exponential scale is linear, so the difference $\Delta(\ln (z))$ or $d(\ln (z))$ is constant, and equal to $\ln (z(1))da$. To better see this, let's derive directly for this example, using $n_v$ for "number of voices":
$$
z(a) = 2^{a/n_v} \rightarrow dz(a) = z \frac{\ln(2)}{n_v} da = \ln(a_0)da \tag{5}
$$
When discretizing, $da=1$, and we end up with
$$
\boxed{f[] \approx \frac{\ln (a_0)}{2 C_\psi}  \Re e \left\{\sum_{k=1}^{K} \frac{1}{a_0^k} W_f (a_0^k, .) \right\} } \tag{6}
$$
same as Mallat Eq 4.67 (note his derivation is for double-integral, but same scaling argument applies here).


Are there accepted ranges of "high", "medium" and "low" visual spatial frequency? I have an image and have been able to quantify its power across spatial frequencies, using Python. What I am interested in doing is creating some measure of how much of the image's information is low, medium and high SF.

I figured a good way to do this would be to take the average power within three bins of frequencies. This leads to two questions:

Is that actually a good way to do it?
Are there agreed upon ranges of "low", "medium" and "high" spatial frequency?

Thanks!

The short answer is yes(1), no(2).

The quantification of the power spectrum is performed via the Fourier Transform over an image. An image is a matrix that describes the distribution of reflectance (or radiation more generally) across the field of view of the camera. It therefore has only a remote "relationship" with what is actually depicted in the image (the physical scene). Image resolution and spatial resolution are two different things. You can use a 50MP camera at the end of a telescope or at the end of a microscope. Both images will have the same image resolution but they will be "looking at" different spatial resolutions (actual space being sampled).

In other words, small objects correspond to high frequencies and large objects correspond to low frequencies but what is large and small depends on the context of the image (and the characteristics of the imaging equipment to an extent). Consequently, what is "high", "medium" and "low" spatial frequency can only be expressed in relative terms and either with a lot of assumptions or constraints.

The approach of dividing the available bandwidth into three bands and estimating the average power in each is, of course, reasonable.

An improvement would be to take a set of representative images, average their spectra and look at the actual bandwidth they occupy and its shape with a reference to the actual objects they represent to then decide at a different possible assignment between the low, med, high characterisations. But this implies specific scenes with a specific field of view, orientation, resolution and any other factor that affects the recorded spatial frequencies.

In case it is of any additional help, there is a concept in image histograms called "Adaptive Binning" where, the width of a given histogram bin is not fixed. The width is decided after an optimisation process over a criterion and it is image specific. For example, one criterion could be that each bin should represent at least N pixels. Although originally not defined over spectra, this adaptive concept might be useful here as well but applied over FFT bins. For more information please see this link.

Hope this helps.


Edge detection in spatial or in frequency domain? I want to ask why everywhere is used to do edge detection in spatial domain and not in frequency domain with FFT of image and apply HP filter. Thank you.

What you claim is simply not true. One step in edge detection is often a simple high-pass (or band-pass) filter,  and those are commonly applied using fast convolution, which is a frequency-domain algorithm.


What approximation techniques exist for the square super-root function? I need to implement an approximation to the inverse of $x^x$, i.e. the square super-root (ssrt) function. For example, $\mathrm{ssrt}(2) \approx 1.56$ means that $1.56^{1.56} \approx 2$.  I'm not as interested in any particular accuracy/bit-depth as I am in understanding what my options are in contrast to more straightforward approaches using power series.

Wolfram Alpha gives a nice symbolic solution in terms of the Lambert W function (i.e. $\ln(x)/W(\ln(x))$ ).  Wikipedia gives the same formula, as well as the equivalent $e^{W(\ln(x))}$. Given that there's a reasonable amount of information on computing $W(x)$ [1] [2], technically that's everything needed to implement something for a variety of requirements.  I know of at least two books that go into extensive detail about approximating $\ln(x)$ [3] [4], so there's even plenty of room to optimize from that direction.

However, I have two questions:

Have approximation techniques specific to this function been published anywhere?
Does it go by another name besides "square super-root" that would make searching for references at a little easier?

Wikipedia/Google has turned up some references dedicated to more general "tetration" functions that include $\mathrm{ssrt}(x)$ as a special case, but most of them seem to be more geared to exploring/defining the general cases.

--

Corless, R.; Gonnet, G.; Hare, D.; Jeffrey, D.; Knuth, Donald (1996), "On the Lambert W function" http://www.apmaths.uwo.ca/~djeffrey/Offprints/W-adv-cm.pdf
Digital Library of Mathematical Functions. http://dlmf.nist.gov/4.13
Crenshaw, Jack W. (2000), Math Toolkit for Real-Time Programming.
Hart, John F. (1978), Computer Approximations.
Chapeau-Blondeau, F. and Monir, A. (2002). Numerical evaluation of the Lambert W
function and application to generation of generalized Gaussian noise with exponent 1/2. IEEE
Transactions on Signal Processing 50, 2160-2165. http://www.istia.univ-angers.fr/~chapeau/papers/lambertw.pdf
Minero, Paul. Fast Approximate Lambert W. http://www.machinedlearnings.com/2011/07/fast-approximate-lambert-w.html

--

Update

After doing some more research over the past few days, I still haven't found the kind of hands-on "Crenshaw style" $[3]$ treatment of $\mathrm{ssrt}(x)$ I was hoping for, but I did find a new reference worth documenting here.  On page three in $[5]$, there's a section titled "Fast Approximation" that goes into great detail about approximating $W(x)$ in the context of noise generation. As an interesting aside, the probability density of "Gaussian noise with exponent 1/2" [in the paper] looks strikingly similar to the histogram in Kellenjb's answer to this question about detecting signal clipping.

In addition, the link given by rwong in the comments $[6]$ is a great resource for actually implementing $W(x)$, and it even links to the author's BSD licensed project called fastapprox, which includes the implementation described.

Some numerical stabs in the dark yielded the following for an iterative approach:

We're looking for the solution y = f(x) where y^y = x.

$$y \ln y = \ln x$$

$$y = g(x,y) = e^{\frac{\ln x}{y}}$$

The value for $y$ is a fixed point of the above equation, and empirically it seems to converge for some values of $x$, but for larger values of $x$ it oscillates or diverges.

Then I tried an approach similar to Newton's iterative square-root:

$$y = \frac{y_{previous} + y_{*}}{2} = \frac{y + e^{\frac{\ln x}{y}}}{2}$$

where y* is supposed to represent a nonconverging but optimistic answer that maintains accuracy if you happen to guess an accurate initial value (in square root y2 = x, it's y* = x/y).

This appears to converge, but very slowly at the low end of $x$ (near $x_{min} = (\frac{1}{e})^{\frac{1}{e}}$ )

It also looks like a good initial guess is $y_0 = \ln(x)+1$.

So I figured maybe there's a better-converging solution:

$y = (1-a) \times y + a \times g(x,y)$ for some value of $a$ that is a function of $x$.

Then I found something interesting.

If I get a converging answer $y$ from the above approach for $y^y = x$, and then compute $y_2 = g(x,y+ \epsilon) = e^{\frac{\ln(x)}{y+\epsilon}}$, it appears as though $y_2-y$ = approximately $\epsilon \times (-\ln(y))$.... e.g. if we had a guess $y_1 = y + \epsilon$ for some unknown $\epsilon$, and computed $y_2 = g(x,y_1)$, then $(y_2 - y) \approx \epsilon \times (-\ln(y)) = (y_1 - y) \times (-\ln(y))$. (Just to clarify, I have no analysis to verify this, but the numbers just popped out of some numerical evaluation I performed.)

Solve for the linear terms in $y$, and you get $y = \dfrac{y_2 + \ln(y) \times y_1}{1 + \ln(y)}$... use $\ln(y_1)$ in place of $\ln(y)$ and you get this iterative approximation:

$$\begin{align*}
    y[n+1] &= \frac{g(x,y[n]) + \ln(y[n])\times y[n]}{1+\ln(y[n])}\\
           &= \frac{e^{\frac{\ln(x)}{y[n]}} + \ln(y[n])\times y[n]}{1+\ln(y[n])}
\end{align*}$$

This appears to work very well, with the initial guess $y = 1+\ln(x)$, and appears to converge within 4 or 5 iterations.

(Someone could probably show that this is equivalent to Newton-Raphson in some way, but I think it's beyond my capability.)


Dirichlet Conditions What is meaning of a signal having a "finite number of maxima and minima during any single period of time"?

The Dirichlet conditions provide sufficient conditions under which a real periodic function can be assimilated to its Fourier series at each point of continuity.

One of these conditions is that the functions have a finite number of extrema  (maxima or minima) in every bounded interval. This means that, in some sense, the function does not oscillate infinitely on some time segment, even with very tiny oscillations. An example of such an oscillating continuous function can be given by:

$$ f(t)=t^2 \cos (1/t)$$

with $f(0)=0$ by continuity. The ripples around $0$ get narrower and more numerous, and are infinitely many, in any non empty interval containing $0$.

The function is evidently continuous, even differentiable. On each interval defined by $$\left[\frac{1}{(k+1)\pi+\pi/2},\frac{1}{k\pi+\pi/2}\right]$$
the function $f$ vanishes at both ends, and since not identically $0$, it attains at least one extremum on this interval. Hence the function has an infinity of extrema on each interval $[0,a]$, $a>0$. The exact locations of the extremas can be computed with the derivative of $f$.


What does one need to know in order to build a software synthesizer? For someone who has some experience in software development, loves his synths and knows the [very] basics of digital audio representation, what are the first steps he should follow in order to program a music synthesizer, DSP-wise?

Get Real Sound Synthesis for Interactive Applications, and download the Synthesis Toolkit.  That will give you practical introductions and quality source code for most digital synthesis methods in use.  Then pick whichever method you find most interesting, and start modifying the corresponding STK code to do what you want.  Then try building something similar from scratch.

In addition, like hotpaw2 suggested, anything from JOS is also great introductory material.  In particular, Introduction to Digital Filters with Audio Applications and Mathematics of the Discrete Fourier Transform both start from the basics, and don't assume a lot of background.

From a digital filter/effects perspective, another good starting point is to check out Faust.  You don't even have to learn their block-diagram signal processing language, you can just use the examples to dump out VST, LADSPA, etc. template code as quick way of getting something running [which you can then modify].

That being said, it's even easier to start out with sound rendering, which is basically the same thing, but without worrying about real-time operation or interactivity.  You just save wav files to disk.  It's a great way to experiment with new techniques before going through the trouble of making them run quickly.


Why OFDM is more appropriate with MIMO than other multiple access method? As read from MIMO-OFDM wikipidia page says

  Research conducted during the mid-1990s showed that while MIMO can be used with other popular air interfaces such as time-division multiple access (TDMA) and code-division multiple access (CDMA), the combination of MIMO and OFDM is most practical at higher data rates.

Why does OFDM has more higher data rates and they combined more advantaged? Thank you

One possible answer is given at the beginning of the text:

  Raleigh suggested and later proved that the processing required by
  MIMO at higher speeds would be most manageable using OFDM modulation,
  because OFDM converts a high-speed data channel into a number of
  parallel lower-speed channels.

Also:

  MIMO-OFDM is a particularly powerful combination because MIMO does not
  attempt to mitigate multipath propagation and OFDM avoids the need for
  signal equalization.


Transmitting over Unknown Fading Channel I have a system with no diversity and the channel's coefficients change every $T = 10\text{ bits}$ transmitted. The channel estimation system has a flaw and detects the channel's coefficients equal to $h$ while their real values are $h+c$ where $c$ is a constant equal to $1$.

Is there any method I can apply to achieve a better BER than simply using the faulty $h$?

I am using BPSK.

The strategy depends heavily your "realistic scenario", e.g. which estimator, which equalizer, which estimation error, etc.

In general scenario, you cannot trust the estimate of $h$, it means that you cannot use coheent detection. Thus try non-coherent detection (detection by energy) : to send bit $0$, use 2 channel uses (2 symbols) $x_0 = [x[0], x[1]] = [1, 0]$ and to send bit $1$, use  $x_1 = [x[0], x[1]] = [0, 1]$. The receive symbols are $y = [y[0], y[1]]$

In the well-known Rayleigh flat fading channel, i.e. $y[m] = h[m]x[m]+w[m]$, the optimal detection rule will be
\begin{align}
\hat{x} &= x_0 \textrm{ if } |y[0]|^2 > |y[1]|^2\\
\hat{x} &= x_1 \textrm{ otherwise}
\end{align}


How do I go from LPC coefficients to a filter polynomial? I've recently begun experimenting with LPC, and while I understand that it works, I'm still slightly lost on why it works.

Specifically, I understand that LPC involves finding coefficients $a_1, a_2, \ldots a_p$ such that for a signal $h$,

\begin{align}
h(n) &\approx \sum_{i=1}^{p} a_i h(n-i) \\ &= a_1 h(n-1)+a_2 h(n-2)+\ldots+a_p h(n-p)\end{align}

In other words, we want to find $a_1, a_2, \ldots a_p$ which will let us linearly approximate the next sample given the previous $p$.

However, when LPC is actually used, these coefficients $a_1, a_2, \ldots a_p$ are usually treated as the coefficients of a polynomial:

$$x^p + \sum_{i=1}^{p} a_i x^{p-i} = x^p + a_1 x^{p-1} + a_2 x^{p-2} + \ldots + a_p$$

Matlab's LPC documentation refers to this as "the prediction filter polynomial", and it's very useful for e.g. finding the formants of a signal (which are simply the zeroes of this polynomial).

But I don't understand the connection between these two uses. Why is it useful to turn the linear prediction coefficients into a polynomial like this? And why does it work?

The prediction polynomial is the representation in the $z$-domain of the first equation you wrote.

When performing LPC, one assumes that the filter consists of only poles and the transfer function is $\frac{1}{A(z)}$ where

$$A(z)=1 - \sum_{i=1}^{M} a_i z^{-i}$$

That's the polynomial you are talking about. Note that it should have only negative powers of $z$, not positive.

As you wrote in your question, one estimates the signal $h(n)$ based on past values of that same signal. The $z^{-i}$ factors in the polynomial represent those delays in the $z$-domain, while the coefficients $a_i$ are just constants and, due to the $z$-transform linearity, they just stay as they are.


Impulse response of forward difference cascaded with one sample delay Below is the excerpt from Discrete Time Signal Processing by Alan Oppenheim.

I don't get how $(\delta[n+1] - \delta[n]) * \delta[n-1])$ becomes $\delta[n] - \delta[n-1]$. The convolution sum operator "*" here is defined as below.
$$y[n] = \sum_{k=-\infty}^{\infty} x[k]h[n-k] $$
How do I substitute the definitions of forward difference and one sample delay in this and arrive at the impulse response the author has given?

For any sequence $x[n]$ you have
$$x[n]\star\delta[n-n_0]=x[n-n_0]\tag{1}$$
because
$$x[n]\star\delta[n-n_0]=\sum_{k=-\infty}^{\infty}x[k]\delta[n-n_0-k]\tag{2}$$
Note that $\delta[n-n_0-k]$ is only non-zero for $k=n-n_0$, hence the result $(1)$.
Now use $x[n]=\delta[n+1]-\delta[n]$ and $n_0=1$.


Help understanding impulse response I am having trouble understanding the impulse response function.
 I have a beam of 10 m long. I place 2 sensors, one at point A and one at point B. I then excite the beam at point E (shown in picture) by tapping it with a hammer. I then repeat this on several different beams of different lengths, dimensions, materials and more. I want to tap the beam to see how the beam responds, using accelerometers at locations a and b, and see if I can estimate things about the beam (example material).

In the case of calculating impulse response, what exactly is my input? I assume it is the excitation signal, however I do not measure the excitation signal until it gets to point A (there is no sensor at point E). I therefore do not know the excitation at the source, it just be with a hammer hitting the beam. The matlab documentation uses inputs as an example, but this does not mean anything to me as I do not know my input.

Also, what does the impulse response actually tell me? How is this useful?

I am using matlagb.
Gracias todos,
Ben

Firstly, for a response to be classified as the response to an impulse. The input first must qualify as an impulse to the system. For ex: a quick clap could be the impulse input to measure the sound impulse response of the room.

In your case, hitting the hammer at point E must be a "good enough" impulse input for all the materials and lengths.

The issue however is complicated in this case as locations A and B are simply outputs when viewed inline with systems theory for the rope system you constructed. Therefore, the excitation at input should be an "effective impulse" at both locations A and B for you to effectively denote the response of both sensors at A and B to be the response of an impulse.

So the essential question here is to first classify what qualifies as an impulse for rope at E and the suseuqent uniformity of excitation across the length.

If the excitation  at E is not uniform at A and B then you would have to define two effective impulses for your rope system. One for the response at A and other at B.

Once you know the impulse response. Then any input can be written as sum of weighted impulses and if the system is linear and time invariant then the output is simply the sum of delayed responses as below for a discrete time signal

$$ x(n) = \sum_{k=-\infty}^{k=\infty} x(k)\delta(n - k)$$

Let the impulse response to this $\delta(n)$ input be $h(n)$ then the output is simply,

$$y(n) = \sum_{k=-\infty}^{k=\infty} x(k)h(n-k)$$

This is possible both because the system is linear as well as time invariant, so if the impulse response $\delta(n)$ is $ h(n)$ then the response for $\delta(n-k)$ is $h(n-k)$