How to choose FFT depth for ADC performance analysis (SINAD, ENOB) I am trying to simulate a model of an ADC and determine its performance.

One of the interesting properties is the ENOB (Effective Number Of Bits), which can be calculated from SINAD (SIgnal-to-Noise And Distortion ratio).

On that SINAD Wikipedia page, there is a PDF document which suggests that the second definition on the SINAD Wikipedia page is the one to use and which I interpreted to be

$$
SINAD = 10 \log_{10} \left( \frac{p_f} {\sum_i{(p_i)} - p_0 - p_f} \right)dB
$$

Where $p_x$ is the power the FFT bin at frequency $x$, $p_f$ is the power of the frequency bin containing the signal frequency $f$ and $p_0$ is the DC component. I calculate the power of each bin by squaring the normalized amplitude. Also note that the sum in that equation runs over all frequencies from $0$ to the Nyquist bandwidth $f_s/2$ with $f_s$ being the sampling frequency of the ADC.

That same document also defines
$$
ENOB = \frac{SINAD-1.76dB}{6.02}
$$

Using my SINAD and this ENOB definition, I compared the output of an ideal $N$ bit ADC model, which matched for my FFT depth.

In another (related) PDF document, it states that the FFT depth must be large enough in order to distinguish between the FFT noise floor and the ADC noise. To my surprise, it also says that the FFT noise floor is
$$
10 \log_{10} \left( \frac{M}{2} \right)dB
$$
below the theoretical signal-to-quantization-noise $SNR=6.02N+1.76dB$, which makes me wonder how $M$ can thus ever be "too low"?

While playing with different ADC models, I have seen that the ENOB does vary significantly when I set $M$ too low, so I am wondering if there is any guide on how to choose $M$ for a desired ADC bitwidth of $N$?

tl;dr: ADC noise is both distributed (approximated white noise) and spurious. A N-point FFT is a bank of N filters each approximating a Sinc function in frequency with an equivalent BW of 1 bin. The power output of each FFT bin is proportional to the total power across that filter's BW. For the case of white noise, as the BW gets bigger, the power out in each bin goes up--- this is the FFT noise floor. Spurs, unlike white noise, are not a distribution and therefore do not go up and down as we change the BW of the filter (assuming the filter is centered on the spur). In order to have enough dynamic range in the FFT to see the spurs, we need to ensure the filter BW is narrow enough so that the total integrated white noise power in the filter is less than the total power from the spur.

More Details:

Quantization noise can be well approximated (modeled) as uniform white noise, meaning the distribution of the amplitude is uniform (with an uncorrelated input waveform and sufficient number of bits, it is equally likely that "truth" is anywhere between the quantization levels at each sampling interval), and the distribution of the energy in frequency is also uniform (white) over the sampled spectrum from $-f_s/2$ to $+f_s/2$.  In practical implementation due to correlation with the test signal, and non-linearities in the ADC, will result in individual frequency tones in the ADC output called spurs. Thus both Spurious Free Dynamic Range (SFDR) and Signal-to Noise and Distortion (SINAD) are specified (among other ADC performance metrics).

The SFDR is the rms level of the strongest sinusoidal input that maximizes the SFDR (typically this level is around 1 dB below where the actual full-scale is) relative to the level of the maximum spur rms level.

Signal to Noise and Distortion (SINAD) is typically used to denote the rms level of the Signal relative to ALL noise contributions (meaning all the spurs, all the quantization noise, etc), again at a Signal level for a sinusoidal test tone that would maximize this metric. DC offset is usually not included.  Since SINAD includes all spurs, including the strongest one, SINAD will always be less than SFDR; however SINAD represents ALL the noise power in the digital spectrum, where SFDR represents a tone at one specific frequency.

After I will show how the formula for SINAD =  6dB/bit + 1.76 dB is derived which will give further insight, but first to answer your question: A key point is that the noise component of SINAD is both at individual frequency locations as spurs and also as distributed noise. Together they are the total noise in SINAD.  An FFT is a bank of filters, each filter separated by $F_s/N$ where $F_s$ is the sampling rate and N is the number of bins. (This is very clear if you do a sliding FFT across a waveform, recomputing a new FFT output after each sample shift such that a waveform is produced at each bin). Further, without modifying the time domain data with windowing or zero padding, each bin has a bandwidth of 1/T where T is the duration in seconds of the time domain block being FFT'd. So there is a filter for each bin and each filter is equivalently 1 bin wide. Also note that as the number of samples in the FFT increases, the shape of each filter in frequency approaches a Sinc function (which isn't a great filter, meaning each FFT output will be sensitive to frequencies in other bins).

This is demonstrated with the simple case of a 4 point FFT in the figure below, where we see 4 filter responses, one for each of the 4 FFT outputs. The output of each "filter" will be proportional to the total integrated power across the spectrum after being filtered in each bin; thus in this example we see how a single tone shows up in all the bins, based on each bin's respective attenuation at the particular frequency the test tone is at.  However each bin has an equivalent noise band width of 1 bin which is equal to $1/T$ as I described earlier. This means that if the input was white noise, an equivalent brickwall filter with bandwidth $1/T$ would have the same total noise power at the output of the filter.  So for the case below, with a white noise input, each of the four FFT outputs would have a noise power term that is 1/4 of the total noise. (If we did successive 4 point FFT's from a larger data set, either block by block or a "moving" FFT where we slide over the sample set one sample at a time, it would be very clear by observing the outputs that their rms level would each be 1/4 of the input rms level, assuming a white noise input).

Every time we double the number of points in the FFT, each filter reduces in BW by half, so for white or similarly distributed noise, the power would reduce by $10log(2)$ = 3 dB.  The power level for a spurious signal which exists at one frequency (and is not a power distribution) would not change. Therefore in order to see a spurious signal, we need to ensure we have enough resolution (narrow bandwidth filters) so that the total distributed noise power in the bandwidth of each bin is less than the noise power of the spur.

Why 6 dB/bit + 1.76 dB?

In the plot below I show the SNR for a "perfect" A/D Converter with a noise distribution due to quantization noise modeled as a uniform distribution (noise that is uniform in amplitude over a height of one least-significant bit). Since this noise is due to quantization noise alone, it is also called SQNR (Signal to Quanitization Noise Ratio).

The first plot shows the distribution (histogram) for a full scale Sine Wave, along with the variance for such a distribution. This is compared to the variance for the uniformly distributed quantization noise. As the next slide slows, it is the ratio of these two distributions (Signal - Quantization Noise) that results in the handy relationship between number of bits and SQNR, importantly for a full scale sine wave:

SQNR = 6 dB / bit + 1.76 dB.

I had actually completed and experiment to see how well this approximation would hold up for lower number of bits with the results graphed below, showing little error with the predicted result over a wide range of number of ibts (7 bits results in 1/11.53 instead of 1/12 which would result in an error in the calculated SNQR of only 0.17 dB.)

Further it is interesting to see that with a large number of bits, the histogram is uniformly distributed as expected, and the quantization noise is fairly "scrambled" relative to the input signal and from sample to sample; resulting in a reasonably white spectral density in frequency (unfortunately frequency plot is not shown).

In comparison, here is a plot at the other extreme with a 1.585 bit converter. Notice that the histogram approaches that of the input signal (the sinusoid), explaining why the early plot showing the variance for quantization noise at different number of bits approaches 1/8 as the number of bits is reduces (approaching that of the input signal!).  Also note the pattern of the quantization noise itself plotted below the signal (showing the error from "truth" to the levels available). Here we see clear correlation both with the input signal and from sample to sample; this means very strong spurs will be present!

The above cases were for a modeled "perfect" A/D converter; more specifically a quantized digital signal. In actual ADC implementation, there will be other analog distortion sources and non-linearities that will also contribute to noise. This gives rise to Effective Number of Bits (ENOB). ENOB converts the actual A/D Converter into an equivalent "perfect" A/D converter: specifically an equality is made between the actual measured or expected Signal to Noise ratio from all noise sources (SINAD) and SNQR which implies quantization is the only noise source, and by solving for the number of bits in SNQRthat would be equal to SINAD  (from 6dB/bit + 1.76 dB as the OP had done) , we determine the equivalent "perfect" A/D converter that would result in the same total noise level. Because the noise is both distributed and spurious, it is important to consider both ENOB (or SINAD, they are both derived from the same result) and SFDR.

Here is also a good reference from Analog Devices that gives further details:

https://www.researchgate.net/file.PostFileLoader.html?id=58871498404854e3dc2f909f&assetKey=AS%3A453975662174208%401485247640061


Continuous and discrete-time raised cosine filter properties in frequency domain In communication systems, the raised cosine (RC) filter is split into root-raised cosine (RRC) filter at the transmitter and the receiver. The combined response of both RRC filters is the RC filter.

If the transmit RRC filter is analog (continuous-time or CT), how the combined response of both filters change if the receive RRC is placed (a) before digital sampling and (b) after?

For example, let $h_{\rm T_{RRC}}(t)$, $h_{\rm R_{RRC}}(t)$ and $h_{\rm RC}(t)$ be the impulse responses of the transmit RRC, receive RRC and RC filters, respectively. Let their frequency responses be $H_{\rm T_{RRC}}(f)$, $H_{\rm R_{RRC}}(f)$, and $H_{\rm RC}(f)$, respectively. Let the symbol period be $T_\rm c$.

Case 1: Receive RRC filter is continuous-time

Here the system is as follows:

Transmit RRC $\longrightarrow$ Receive RRC $\longrightarrow$ Sampling at rate $1/T_c$

The combined response of the filters before sampling is
$$H_{\rm T_{RRC}}(f)H_{\rm R_{RRC}}(f) = H_{\rm RC}(f) \tag{1}$$

After sampling, the RC filter has the following desired property in time-domain
$$h_{\rm RC}[n] = h_{\rm RC}(nT_\rm c) = \begin{cases} 1, \: n = 0 \\ 0, \: \text{otherwise}\end{cases}\tag{2}$$
$$\textrm{or,}\quad h_{\rm RC}(t)\cdot\sum\limits_{k=-\infty}^{\infty}\delta(t-kT_\rm c) = \delta(t)\tag{3}$$

In frequency domain, this becomes

$$H_{\rm RC}(f)\star\frac{1}{T_\rm c}\sum\limits_{k=-\infty}^{+\infty}\delta\left(f - \frac{k}{T_\rm c}\right) = 1\tag{4}$$
$$\textrm{or,}\quad\frac{1}{T_\rm c}\sum\limits_{k=-\infty}^{+\infty}H_{\rm RC}\left(f - \frac{k}{T_\rm c}\right) = 1.\tag{5}$$

Case 2: Receive RRC filter is discrete-time

Here the system is as follows:

Transmit RRC $\longrightarrow$ Sampling at rate $1/T_\rm c$ $\longrightarrow$ Receive RRC

Here, the sampled transmit RRC response is
$$\frac{1}{T_\rm c}\sum\limits_{k=-\infty}^{+\infty}H_{\rm T_{RRC}}\left(f - \frac{k}{T_\rm c}\right).$$

Would the following hold true now:
$$H_{\rm R_{RRC}}(f)\cdot\frac{1}{T_\rm c}\sum\limits_{k=-\infty}^{+\infty}H_{\rm T_{RRC}}\left(f - \frac{k}{T_\rm c}\right) = 1?\tag{6}$$

-ryan

Root Raised Cosine has a bandwidth $B=\frac{1+\alpha}{2T_c}$ which is $(1+\alpha)$ times more than the minimum required bandwidth set by Nyquist no-ISI criterion. To satisfy the Nyquist sampling theorem now, one must sample at a rate of $2B$. In time domain, we can say that we have to have more samples/symbol than 1. In practice, and to facilitate the synchronization blocks down the chain, this sampling rate is selected such that we have an integer number of samples/symbol. The smallest integer greater than 1 is 2. In your eq (6), you are sampling at a rate of 1 sample/symbol which induces aliasing. Hence, the results of eq (5) and eq (6) will not be equivalent. For 2 or more samples/symbol for incoming signal and assuming no analog distortions, eq (5) and eq (6) will become equivalent after downsampling to symbol rate.

The best way to understand is to imagine a figure of RRC in frequency domain with aliases at integer multiples of symbol rate $1/T_c$.


Can a software defined radio with a frequency range of 300 MHz – 3.8 GHz detect AM (535 to 1705 KHz) and FM signals (88 to 108 MHz)? I am interested in knowing whether BladeRF SDR with a frequency range of (300 MHz - 3.8 GHz) detect AM and FM frequencies? Also, how important is the sampling rate of the SDR. For example, BladeRF has a sampling rate of 80 KSPS - 40 MSPS. What does this number signify? will every signal going out/in the SDR be sampled with this sample rate? Similarly, what does the maximum bandwidth signify?

For those given specs, without additional hardware (upconverters) or modifications to the SDR circuit, it's not possible to detect AM or FM (88-108 MHz) range of frequencies. As you might have also seen, there are frequency range extenders (especially on the low frequency range) available as kits.

Also, the sampling rate determines the maximum bandwidth of simultaneous detection (or spectrum analysis if you prefer to use the SDR for that purpose). So, the bigger the sampling rate the larger will be the spectral bandwidth that you can capture instantaneously. This means, more radio channels are observed at once and more carriers are searched without changing the central tuning frequency of the SDR tuner chip. Which makes it faster in spectrum scanning applications.

However, using a large bandwidth (high sampling rate) mode of SDR, might bring the problem of transmitting the data from SDR to PC, especially on a USB 2.0 interface. (this is not a problem on most ubiquitous SDR receivers used at 2.4 MSPS rates) Also more processing and interrupting than necessary (provided a narrow bandwidth radio channel is to be demodulated) is done at those high sampling rates, to decode just a small portion of the analysis bandwidth.


Direct Form 2 IIR filter input values In the filter graph below:

where circle with 'X' inside is a multiplication
circle with 'Σ' inside addition
and the triangle is a delay element

What could be the values x(k) in this type of filter? There is a rule? There is a range? I can understand the matrix of coefficients but i don't know about the x(k) values.

$x[k]$ is the input signal to the filter; it could be anything. There is no "rule" for what signal you can input to a linear filter. With that said, the stability of filters is usually talked about in a bounded-in, bounded-out (BIBO) sense. If the input signal is unbounded, then it's likely that the output signal will be also. So if anything, the "rule" that you're looking for would be to only apply signals $x[k]$ whose magnitude $|x[k]|$ is bounded below some value $B < \infty$.


How does shifting a signal in the time domain affect its frequency domain? Suppose the signal is shifted by dt (signal 'starts' later, say after 1s instead of 0s), does that correspond to a positive or a negative phase shift df in the frequency domain?

There are certainly very detailed answers to this but I am still having trouble, especially with the sign conventions.

Thanks

According to the Fourier's transform properties, the Fourier's transform of $f(t+dt)$ would be $F(f)e(j2\pi fdt)$
So, basically, the spectrum (Fourier's transform magnitude, representing the frequency content of your signal) does not change. However, f Fourier's now has a phase-shift proportional to the frequency.

To answer your question : the phase-shift's sign is the time delay's sign.


Discrete Fourier Transform in Signal Processing - Interpreting graphs of transformed signals Given above are the real parts of the signals I to IV. Which of the following statements are correct?

(i): Signal III is the result of the discrete Fourier transform of signal I. The associated imaginary part is 0.

(ii): Signal III is the frequency spectrum of an oscillation with a frequency.

(iii): Signal IV is the frequency spectrum of an oscillation with 2 different frequencies.

(iv): Signal IV is the result of the discrete Fourier transform of Signal III. The associated imaginary part is 0.

(iii) seems be to correct. Spectrum of two frequencies.

The other choices can be verified,

(i) Fourier transform of $f(t=0)=1$ is 1. So it cannot be graph (iii)

(ii) Signal (iii) is sine wave and there is no spectrum associated with that.

(iv) It's Fourier transform of sine wave (iii) which is

$$f(t)=0.5\cdot \sin(2\pi Ft)$$

and it's transform is

$$\frac{\pi}{i}\left[\delta\left(\omega-2\pi F\right)-\delta\left(\omega+2\pi F\right)\right]$$

Imaginary part is not zero


Theoretical question on Fast Fourier I have a question regarding Discrete Fourier Transform and Discrete Time Fourier Transform.

Suppose that we have a sequence $x(n)= \{1,2,3,3,2,1\}$.

So is it possible to reconstruct the Discrete Time Fourier Transform $X(e^{j\omega})$ from the Discrete Fourier Transform $X(k)$ ? And how?

In my opinion , it is possible, although I am not sure if something like this would be any good. But in order to reconstruct it we probably need as much Discrete Fourier Transform samples as possible, right?

What do you people think on this matter?

This is (trivially) possible because all Fourier transform like operations are ideally invertible. So what you ask for is actually easily realized by:

$X(e^{j\Omega}) = DTFT\{ IDFT\{X(k)\} \}$

In formulae this is

$X(e^{j\Omega}) = \sum_{l=0}^{N-1} \sum_{k=0}^{N-1} X(k)e^{2 \pi j kl/N} e^{-j\Omega l}$

which can be a bit simplified to

$X(e^{j\Omega}) = \sum_{k=0}^{N-1} X(k) \sum_{l=0}^{N-1} e^{jl (2 \pi k/N-\Omega)}$

or, finally

$X(e^{j\Omega}) = \sum_{k=0}^{N-1} X(k) \frac{e^{2 \pi j k- jN\Omega}-1}{e^{2 \pi j k/N-\Omega}-1}$.

EDIT: The question on whether you need all samples X(k) lies in the behavior of the interpolation function $I(x) = \frac{e^{j N x}-1}{e^{j x}-1}$. This is because the above result is then equal to

$X(e^{j\Omega}) = \sum_{k=0}^{N-1} X(k) \, I(2 \pi k/N - \Omega)$.

I(x) has some "sinc like" behavior, so the question on how many X(k) you need, is really when I(x) has decayed to a value "small enough".


How and why does downsampling reduces resolution in time in my Course we talked about filter banks and the prof said that we get reduced resolution in time because of the downsampling after the band-pass filters. Is this because we now have fewer sample points per second to look at the narrowband signals?

Anytime you take away "actual" data samples away by donwsampling, yes, you are loosing on finer granular data in time domain and hence reducing resolution. However, consider, the following case, an input signal is first up sampled by 2(interpolation) and then down sampled by 2, does it decrease resolution of actual data, not quite, because the upsampled new data points were anyways interpolated versions of the original data. So you have in total not lost any resolution of the actual data. So when we talk about resolution we should always be referring to resolution with reference to the actual signal data points

It's always helpful to keep in mind what happened before the down sampling.

Almost all the time in a system, a downsampling would be happening on an "oversampled signal"(for noise reduction purposes) or preceded by upsampling for rate conversion. In both cases, the data points that would be throwing off were in any case "not required" for processing.

In summary, yes, any time you throw away actual samples of the signal you do loose on information in time domain, but as illustrated above its important to keep in mind that what is being down sampled are actual or artificial data.

Side note on Frequency Domain: Any time we down-sample a signal, we are essentially limiting the signal bandwidth that can be represented unambiguously


Recorded an audio file using microphone via Windows Sound recorder, but has 2 channels that are slightly different. Why are they different? I recorded sounds using a plug in microphone on a Windows 7 laptop. I used the Windows Sound recorder to record the sounds. It was only one microphone, but when I load the sound file into Audacity, it loads 2 channels. This is to be expected, but the sound files are slightly different (subtracting them from each other in MATLAB doesn't result in a zeroed signal).

Shouldn't they be copies of one another?  I didn't use any special effects etc.

Thanks

I am not %100 sure, but it sounds like you are recording the file in stereo. If this is the case, then Windows Sound recorder is taking stereo input, which means it's accessing 2 AD converters which have the same input applied to them. The difference between the two signals in this case would be the noise and any offset between the two converters.

I would see if Windows Sound Recorder could make a mono recording, or else use another program to make the recording.


Splitting the Unit Delay Please refer to the paper Splitting the Unit Delay - Tools for fractional delay filter design by Laakso, Valimaki et.al.
I am not able to visualize how the fractional delay is obtained by resampling the shifted version of the impulse response of the ideal filter as described in the following para.

When the desired delay D assumes an integer value, the impulse
response Eq. 12 reduces to a single impulse at n = D, but for
noninteger values of D the impulse response is an infinitely long,
shifted and sampled version of the sinc function (Fig. 3).
Unfortunately, the ideal impulse response is not only infinitely long
but also noncausal, which makes it impossible to implement it in
real-time applications.
Equation 12 gives an answer to the original problem, i.e., where the
delayed signal value should be placed as it cannot be put "between the
samples." In the ideal case, it is to be spread over all the
discrete-time signal values, weighted by appropriate values of the
sinc function.

Also, suppose the minimum unit delay supported by a system is D, can we derive fractional delays lesser than D with this technique or only greater than D? Does this in any way overcome system limitation w.r.t minimum tap delay resolution(splitting the unit delay as the name suggests)? Or is it  only useful for implementing fractional or irrational rate conversions within bounds of the hardware?

For example a delay of 0.75 sampling periods would have an impulse response like this (red squares sampled from the blue delayed sinc):

Time (horizontal axis) unit is sampling periods. That kind of a filter is not causal as there are non-zero samples at negative times. None of the samples is zero, no matter how far left or right you go. That's what they meant by spreading all over.

About your second question, I have no access to the paper but I know that such all-pass fractional delay filters can be made arbitrarily good (approaching sinc), and the latency increases as you increase filter order. If you want very short latency, you then will also have a very low-order filter, which is of very poor quality as it is completely missing the left side of the sinc. It is a tradeoff between quality and latency, and you probably don't want to trade away all quality.


upsampling by huge factor Assume i have a discrete signal, which i want to up-sample by factor 100.
Up-sampling using poly-phase algorithm sounds like a bad idea, in terms of cache low-efficiency due to low space-locality of the data.

Do you have any suggestions?

Thanks,

Itay

Upsampling by 100 can be factored into upsampling by 2X, 2X, 5X and 5X for each output block (output blocks sized to be data cache resident), thus requiring a smaller phase table for each FIR interpolation kernel.

Some output points may either need to be overlapped across blocks or re-calculated, depending on computation speed vs cache load latency.


How do I keep a value between known bounds in an open system? I detect points in images with an algorithm. The only parameter is a threshold T that I can modify at each image. Such images are from a continuous video, so each one is only slightly different from the previous one. I don't know a priori how many points I get for a determined threshold, as it also depends on the content of the processed image. The whole system needs at least 20 points to continue running, but the more points I have, the more complex the process is, therefore the slower the application. I have made a scheme to represent this:

I need to maintain nb over 20 and if possible under 50 to avoid over-computation. I determined an "ideal" value at 40. I begin with T = 30, but it's only empirical. Usually, T = 10 gets too many points and T = 60 too little.

My question: How do I filter T depending on past values nb, in order to keep nb around 40 for each new image?

I precise I'm really not an expert in filtering, as you may have guessed. I'm not even sure it's a filter I need here.

EDIT: as requested in comments, here is a plot of nb for five different thresholds. Content is always the same (fixed camera), therefore variations are mostly due to camera noise.

In your plot of nb vs. image it appears that there is some high-frequency "noise" (no surprise there), and that the windowed mean is surprisingly stable.  For the purposes of this answer I am going to assume that the high-frequency noise is fairly representative of the general case, but that there will also be some slow-to-medium trends as well.  It seems like there would have to be if the images are changing.  The frame rate should generally be faster than the changes, though, so the trends should be relatively slow.

Normally I would suggest a low-frequency filter to get rid of the high-frequency "noise", but in this case I don't think that would be wise.  For any given image the next image is very likely to be highly correlated, so I think it makes sense to adjust somewhat if an image has a large change in nb.  I suggest two different solutions.

1) Use a simple moving average of the past n nb values.  Performance is obviously important in this application, and a moving average is computationally cheap.  A moving average gives some stability, which is important given the high-frequency noise, but also is responsive to changes in the data.  You can trade responsiveness for stability by decreasing or increasing n.

Once you have the nb moving average you use that to produce the "error" which is your target nb minus the moving average nb.  Thus, if the average nb is too low, the error is positive, if the average nb is too high, the error is negative.  You use the error to adjust the threshold.  Typically the new threshold would be

$threshold = threshold - \alpha * error$

$\alpha$ needs to be determined either by analysis of your problem or by empirically testing.  You substract the error result from the threshold because a positive error means you need nb to go up, which means that you need the threshold to go down.

2) Use a PID controller.  The PID controller would produce a threshold value for every image, and use the target nb vs. actual nb error (similar to the moving average error except it's an instantaneous error, not an average) to update its threshold for the next image.

The "integral" portion of the controller will give you stability, while the proportional/derivative part will give you responsiveness.  You would need to make sure that the loop itself is stable, and test the constants on different videos, but this would be a very powerful way of adjusting the threshold.  It also offers the possibility of making your app smarter by changing the loop constants dynamically, if needed.

Solution 1 is easier to implement and is a little cheaper computationally, solution 2 is more powerful but also more complicated.  I would see if the moving average is good enough.  If it is, great, if not, implement the PID controller.

EDIT: Tried to make the two solutions a little clearer.


Applications for the Maximum Phase Filter A maximum phase digital filter has all zeros outside the unit circle, and has the maximum phase and therefore longest delay for a given magnitude response.

Besides the possibility of needing a filter with the longest delay possible, are there other practical applications and uses of the maximum phase filter?

I understand that we can deconvolve linear-phase filters into a minimum-phase filter and its reverse a maximum phase filter, but looking for further utility and interesting use cases for maximum-phase filters specifically.

we can deconvolve linear-phase filters into a minimum-phase filter and its reverse a maximum phase filter

We can formulate this more broadly. Any LTI system can be split into a cascade of it's minimum phase filter and an all pass (which is indeed a maximum phase filter). So,
$$H(z) = H_{min}(z) \cdot A(z)$$

where $H_m(z)$ is the minimum phase filter that matches the magnitude $|H(z)|$ and $A(z)$ is an allpass , i.e.

$$A(z) = k \cdot \frac{\prod_{n=0}^{N-1}z-1/p_n^*}{\prod_{n=0}^{N-1}z-p_n}$$

The proof for this is straight forward: Let's assume that H(z) has a zero $q_k$ outdside the unit circle. We can simply factor this out as

$$H(z) = H_k(z) \cdot (z - q_k) = H_k(z) \cdot (z - q_k) \cdot \frac{z - 1/q_k^*}{z - 1/q_k^*}$$

$$ = H_k(z) \cdot (z - 1/q_k^*) \cdot \frac{z - q_k}{z - 1/q_k^*} = H_{k,min}(z) \cdot A_k(z)$$

where $H_{k,min}(z)$ is the version of $H(z)$ where the zero $q_k$ has been made minimum phase by inverting it to $1/q_k^*$  with $|1/q_k^*|  < 1 $ and $A_k(z)$ is a first order allpass with a zero at $q_k$. Note that the magnitude stays the same, i.e. $ |H_{k,min}(z)| = |H(z)| $

That type of splitting an LTI system into it's maximum and minimum phase part is can be quite useful.

Inverting: The minimum phase part is invertible and the inverse is also minimum phase. The inverse of the maximum phase is simply it's own time reverse. That's not directly invertible (since the inverse is non causal), but you easily "eyeball" how much bulk delay you need to add to get enough of the phase response correct.

Audio: Linear phase FIR are mostly useless in audio. Getting decent resolution at low frequencies requires many 1000s of taps and the computational complexity and the latency are often prohibitive. Equalizer or cross over design is often split into two steps: minimum phase filter to get the amplitude right and than allpass filters to do time/phase alignment as needed.

Building blocks: Allpass filters are very useful building blocks to build more complicated structures. For example you an build odd order Butterworth cross overs as sum + difference of two allpasses and as a bonus you get a signal that's phase matched to the perfect reconstruction path. Warped FIR filters are another useful topology where allpasses are quite handy. Beamformers need more or less independent control of amplitude and phase.

So in general in an application where you want to control phase and  amplitude more or less independently, splitting into minimum (for amplitude) and maximum (for phase) can be quite useful.


Is a cubic Lagrange interpolation tensor product the same as bicubic interpolation? I just implemented some interpolated texture sampling by sampling the 4x4 nearest pixels then doing Lagrange interpolation across the x axis to get four values to use Lagrange interpolation on across the y axis.

Is this the same as bicubic interpolation or is it different?

Webgl implementation here:
https://www.shadertoy.com/view/MllSzX

i don't need the points, but the Community keeps resurrecting unanswered questions (or answered questions that have no accepted answer) until i go batty seeing them on the list, time and time-again.

:-\

i think the bicubic interpolation is more likely 3rd-order Hermite polynomial than 3rd-order Lagrange polynomial interpolation. the former guarantees continuity of both the zeroeth and first derivative. the latter only guarantees continuity of the zeroeth derivative (the interpolated function itself).

this previous answer deals pretty completely with the necessary math of it.


DWT edge/indexing problem I'm implementing a very basic forward and inverse DWT, based on Mallat's book. Nothing fancy, simple inefficient matlab code (proof-of-concept more than runtime necessary, so I'm not focusing on optimisation). It almost works. I'm getting a distortion along the top and left edges. I'll post link to (need 10 points to post) some pictures of lichtenstein to exhibit the behaviour.

http://postimg.org/image/igy8cv7hn/
From top-to-bottom:
* After 5 iterations.
* After some unpacking iterations, you can see the distortion along the top and left edges growing.
* Unpacked and distorted.

I'm using symmetric boundary conditions (I've implemented both whole and half, it doesn't make much of a difference) and Haar. Although using CDF has the same effect (the distortion is a little different, but same location along the entire edge). Since the rest of the picture is transformed and inverted perfectly, there is most likely nothing wrong with the implementation as a whole. It must be some kind of index offset. It is clear that the expansion of the distortion is dyadic - it doubles in size for every step of the inverse loop. Two days and lots of index changing has led to no breakthroughs.

I've pasted the code here (I realise it's unlikely people will want to preuse 100 lines of foreign code, but I'm a little desperate). It's pretty simple to run (although like I said it's slow - can take up to a minute to run):

runbackward(runforward(rgb2gray(imread('lichtenstein.png')), 5),5)

Any code-based help or high-level diagnostic ideas would be appreciated.

Bjorne, I am not exactly sure what you are trying to do but I can assure you that the simplest of DWT code is not more than 10 lines in Matlab. Here is how to do it:

approx = zeros(1, cols);
details = zeros(1, cols);

%for Haar Wavelet, Use filterbanks in Matlab, google it!

hi = ...
lo = ...

for i = 1: Level

   a = conv2(signal, lo);
   d = conv2(signal, hi);

   detail() =  downsample(d);
   approx() =  downsample(a);
   signal = downsample(a);

end

That is pretty much everything you need to compute DWT. As for the values between the parenthesis in detail and approximation, just fill in the indexes.


What is the story behind the story about SIFT descriptor? The following is from Lowe 2004 paper ( http://www.cs.ubc.ca/~lowe/papers/ijcv04.pdf ).

One obvious approach would be to sample the local image intensities around the keypoint at the appropriate scale, and to match these using a normalized correlation measure.
However, simple correlation of image patches is highly sensitive to changes that cause misregistration of samples, such as afﬁne or 3D viewpoint change or non-rigid deformations. A
better approach has been demonstrated by Edelman, Intrator, and Poggio (1997). Their proposed representation was based upon a model of biological vision, in particular of complex
neurons in primary visual cortex. These complex neurons respond to a gradient at a particular
orientation and spatial frequency, but the location of the gradient on the retina is allowed to
shift over a small receptive ﬁeld rather than being precisely localized. Edelman et al. hypothesized that the function of these complex neurons was to allow for matching and recognition
of 3D objects from a range of viewpoints.

I am trying to understand SIFT descriptor. I understand the previous stage (keypoint detector).

I don't know why it is implemented that way. I want to know the story behind the story.

The descriptor obtained from a $64\times 64$ neighborhood of interest point at the obtained scale.

It will divide this $64\times 64$ region to $16\times 16$ patches which lead to 16 patches.

For each patch we calculate the gradients and then find the dominant direction of gradients(which has some details), then taking the dominant direction as the reference direction we will divide the 360 degrees to 8 angular region each has 45 degrees, then sum over the magnitude of each gradients which lie in each angular region.

We could consider this as distribution or 8 bin histogram of gradient direction (considering strong gradients has more information we have to use them with higher weight in calculation of distribution so we use their magnitude as their weight which leads to sum over their magnitude). Then we will normalize these histograms.

At the end for each patch we have an 8 bin histogram and we have 16 patches which leads to 128 number descriptor.

By finding dominant direction our descriptor also becomes rotation invariant.
By using gradients our descriptor becomes invariant regarding to baseline illumination and by normalizing obtained histograms our descriptor becomes invariant to the contrast of image.


Convolving two signals I saw a video where this guy used a program to do a frequency analysis on a voice signal and a sawtooth wave (I'm assuming this was FFT). Then he saved the plots as images and combined them pixel by pixel through multiplication using photoshop. He then put this picture back into the program and it did the inverse transform, turning it back into a sound. He said that this was an implementation of a vocoder but the low quality made it hard for me to tell if he was correct.

If so, then couldn't one simply implement a vocoder as a convolution operation? You would just have your two signals then select a window from each of these of the same width at the same position, then perform the convolution operation on these two windows (and probably use a window function as well, such as Hanning). You would, of course, have to do this for every sample, so you would be doing this as many times as you have samples in your tracks (and the windows would sometimes lie partially outside of a track, so they would have to be zero padded).

This seems like it might work because the convolution theorem says convolution in the time domain is itemwise multiplication in the frequency domain, so if it doesn't implement a vocoder, it at least implements exactly what the man in the video was doing (at higher quality). And, I'm not just asking this question blindly, I actually tried it. I get a very cool voice effect but I'm not sure it's the same as a vocoder. In fact, it sounds nothing like a vocoder. If some kind individual could tell me exactly what is happening here, I would be very grateful.

From your description, here is what is happening in the video:

The short-term Fourier transform (aka spectrogram) of the signals is computed. The output of this operation is a matrix of complex values, which cannot be represented as images. Thus, the magnitude or the square of the magnitude is extracted to yield a single positive real value converted into a pixel intensity.
The magnitudes (pixel intensities) are multiplied.
The inverse short-term Fourier transform of the product is synthesized, presumably with made-up phases (or with the original phase of the carrier).

There are so many reasons why this has nothing to do with a convolution. I'm sure other members will point a few more:

What is multiplied are magnitudes or their squares; but not the actual complex values produced by the Fourier analysis.
Even if the quantities that were multiplied were the actual complex values, keep in mind that when dealing with discrete signals of finite length, the inverse discrete Fourier transform of the product of the discrete Fourier transforms is the Circular Convolution, not the convolution.
And still... If you split two signals $x$ and $y$ into blocks of length N, and compute the pair-wise circular-convolution of length N of these blocks, you'll get something very different from $x * y$.

Convolving two audio-signals is a rather meaningless operations. One usually convolve an audio signal with the impulse response of a system; and when convolution is performed in the frequency-domain, the details are trickier than just pair-wise multiplications of FFTs (there's some blocking issues; + necessary overlap add or save).

What you saw would be more accurately described as spectral cross-synthesis. And there is one major reason why this is a very different beast from a vocoder. The goal of a vocoder is to apply the spectral envelope of the modulator onto the carrier. I emphasize on spectral envelope, because when you use a vocoder to make a stack of sawtooth waves say "hello", you apply the formants and overall loudness envelope of the original speech signal to the saw waves, but the last thing you want is the pitch information of the speech signal to get involved. The analysis filter bank of the vocoder should be designed to make abstraction of the individual spectral peaks in the modular signal - what matters is the rough spectral envelope - the location of the bumps (formants). This is why a vocoder doesn't need more than 12-30 channels - too few channels and it doesn't capture the formants, too many channels and it starts capturing pitch-dependent fine spectral peaks. Just like the features used in speech recognition...

Let me give you an example. Let us say you have a modulator speech signal with a 200 Hz f0 ; and formants at 1kHz and 1.5kHz - its spectrum is a sequence of narrow peaks at 200 Hz, 400 Hz, 600Hz, 800 Hz, 1000 Hz, 1200 Hz, 1400 Hz, 1600 Hz, 1800 Hz ; with the peaks at 1000 Hz and 1400 / 1600 Hz emphasized (formants). Let's say your carrier signal is a 140 Hz sawtooth - the spectrum is made of narrow peaks at 140 Hz, 280 Hz, 420 Hz... with a decreasing $1/n$ amplitude. If you go with what you described in your question (product of STFTs), you wouldn't have much left because the two spectra have very little overlap - the first common frequency in the sequence would be at 1400 Hz! What you want to do is somehow capture that the modulator spectrum has a bump at 1kHz and 1.5kHz, and use this information to boost the carrier spectrum in this frequency area. That's how vocoders work.

To do so, and if you really want to go the STFT route, an option would be to smooth the modulator spectrum with a rather wide kernel before doing the multiplication - so that what you are really doing is applying the spectral envelope of one signal onto the other one. This would be akin to applying a motion-blur on the Y axis in Photoshop using this silly image manipulation presentation...


How to perform iris normalization? I have done with iris localization. Now I need to perform iris normalization, but I don't have any idea about it. How to implement iris normalization please?

Please check the following for Iris Normalization:

  Nithyanandam. S, Gayathri K.S., Priyadarshini P.L.K. - A New IRIS Normalization Process For The Recognition System With Cryptographic Techniques

  Arvacheh E. - A Study of Segmentation and Normalization for Iris
  Recognition
  Systems


Are all exponential functions eigensignals of LTI systems? I know that complex exponential functions are eigensignals to LTI systems. Do these include real exponential functions? E.g. $2^t, e^t, ...$
Thanks for the help!

Complex exponentials are eigenfunctions of LTI systems because they are eigenfunctions of the convolution operator:
$$\begin{align}e^{j\omega_0t}\star h(t)&=\int_{-\infty}^{\infty}h(\tau)e^{j\omega_0(t-\tau)}d\tau\\&=e^{j\omega_0t}\int_{-\infty}^{\infty}h(\tau)e^{-j\omega_0\tau}d\tau\\&=e^{j\omega_0t}H(j\omega_0)\end{align}\tag{1}$$
where $h(t)$ is the system's impulse response and $H(j\omega_0)$ is its frequency response (evaluated at $\omega=\omega_0$).
Now try to prove in the same way that the convolution of an impulse response $h(t)$ with a real-valued exponential results in a scaled version of that exponential. Then draw your conclusions.


Diagonal vs non diagonal processing for audio denoising I am going through a technique for audio denoising by block thresholding. Here is the paper:

http://www.di.ens.fr/~mallat/papiers/AudioBlock.pdf

According to the definition in the paper for diagonal estimators,
Simple time-frequency denoising algorithms compute each attenuation factor only from the corresponding noisy coefficient and are thus called diagonal estimators.

and for non diagonal estimators:
To reduce musical noise as well as the estimation risk, several authors have proposed to estimate the a priori SNR with a time-frequency regularization of the a posteriori SNR . The resulting attenuation factors thus depend upon the data values for in a whole neighborhood of (l,k).

I am not quite able to comprehend these definitions, and the difference.
The paper also claims that non-diagonal processing is better than diagonal processing techniques.  Could someone explain the difference in a layman's terms?

I find the diagonal term a bit confusing. Sometimes, for wavelets, people separate scalar and vector or block thresholding. [EDIT] Now I believe I understand it. If you put all your samples or coefficients in a long 1D vector, then multiplying it by a diagonal matrix yield a scalar operation. With non-diagonal, you allow cross-terms.

For structured signals,  or images, or transformed signals (in a wavelet or a time-frequency domain), neighboring samples or coefficients often possess some local correlations, that you want to exploit in denoising. So assume you have an observation $r(k)$ at location $k$ (could be multidimensional), and look for an estimator $\hat{s}(k)$ of some underlying signal $s(k)$, under a noisy model  $r(k) = s(k)+n(k)$. In the original domain, you could for instance apply a scalar mask $f(k)$ on the observation: $\hat{s}(k) = f(k).r(k)$, especially if you know that at some location you have signal, at other noise. Windowing is a typical example. More generally, you can apply a function $f$: $\hat{s}(k) = f(r(k))$. If you assume local correlations, you can borrow information for close samples. For instance, average filtering can be $\hat{s}(k) = \frac{1}{3}(r(k-1)+r(k)+r(k+1))$.

In other words, $f$ now apply over a domain, gathered in a vector $v(k)= [r(k-1),r(k),r(k+1)]$. It can be called Reference Observation Vector (ROV). Then  $\hat{s}(k) = f(v(k))$. So average filters are sort of non-diagonal processing.

Traditionnaly with wavelets or time-frequency, the transform compacts or sparsifies noisy data, and data features tends to become relatively higher than  the noise, compared to the original domain. So the previous idea of a scalar mask works better. You can choose $f$ to keep the value of $r(k)$ if $|r(k)|$ is above a threshold $t$, and set it to zero if below. This $f$ is discontinuous.  let us write it $f(x) = \mathrm{sign}(x) (|x|-t)_+$
This is called scalar wavelet shrinkage or thresholding, scalar because you apply it on a single wavelet coefficient, illustrated in the center plot.

This idea was already used in the Fourier domain for noise reduction in speech. You can have smoother transistion, if you attenuate a bit the higher values, which give you soft thresholding, and you have many more shrinkage functions (garrote, firm) for scalar thresholding.

The next idea is to combine all of the above. You can modify one coefficient with respect to a certain norm of a surrounding $\|v(k)\|$, as:

$\hat{s}(k) = f(\|v(k)\|)r(k).$.

In wavelets, the ROV $v(k)$ can be more complicated. It can gather neighbors of the coefficients, but also coefficients from other scales.

The rationale is that, if a coefficient is small, it might be useful to preserve it anyway. If its ROV has high norm (think of energy), the coefficient might small by mistake, and is not noise.

The nice think is that the Stein Unbiased Risk Estimator theory allows you to estimate something on the signal $s(k)$ that you do not know by using only the observation $r(k)$, using some technical conditions, which allows you to find, automatically, the best threshold $t$, or shape paramaters of the function $f$. Such methods yielded almost the state-of-the-art on wavelet denoising. I am not familiar with audio. However, for multispectral images, here is an example:

a satellite image of the port of Tunis:

its noisy version, with heavy Gaussian noise:

the denoised version:

While not perfect, it was able to recover details, almost invisible from the noisy image, with very few parameters to set (only the shape of the ROV). This answer is a short version of what you can read in  A Nonlinear Stein-Based Estimator for Multichannel Image Denoising, 2008.

Why can  block denoising be better that scalar?  Because scalar is a specific case of block (non-diagonal) processing, if you take the ROV to be restricted to the coefficient $r(k)$ itself. If you choose the time-frequency or time-scale method well, and the other parameters (block shape and size), the block processing is generally (much) better.


Digital Anti Aliasing Filter For Waveform Band Limiting I am working on some waveform generation software and I am trying to figure out what the best type of filter would be to use a an Anti Aliasing filter for my algorithm.
I am generating my waves in a "Raw" mathematical way, meaning that i am creating a ramp for a saw wave and my square wave consists of pure 1's and 0's. The issue is that this technique is inherently not band limited. I would like to implement low pass filter into the algorithm in order to filter of the harmonic being generated that are higher than the nyquist limit.

What would be the ideal filter design for this? What kind of roll of would i need?

I am doing this in c++ but i don't neccesarily need an answer in c++. just an equation or some resources would be good as well.

Thanks

edit:
I am not changing the sampling rate at any point in the chain. But i am still experiencing some aliasing somehow.

here is my code:

double _DSP::Saw_Wave(double* _frequency,_DSP::Saw_Data* _data){
    double _val = _data->_phasor;
    _data->_phasor += 2.0*(1.0/(_DSP::Sample_Rate::Samples_Per_Second()/ *_frequency));
    if (_data->_phasor>_data->_tolerance) {
        _data->_phasor-=2.0;
    }
    return _val ;
}

typedef struct Saw_Data{
        double _phasor = 0.0;
        double _tolerance = 1.0;
    }Saw_Data;

where _DSP::SampleRate::Samples_Per_Second() is the sampling rate.

the code generates the proper wave and frequency but has unwanted low frequency information when played at higher frequencies.

I am generating my waves in a "Raw" mathematical way, meaning that i am creating a ramp for a saw wave and my square wave consists of pure 1's and 0's. The issue is that this technique is inherently not band limited. I would like to implement low pass filter into the algorithm in order to filter of the harmonic being generated that are higher than the nyquist limit.

You can't do this.  You've already aliased the frequencies when you generated the waves in that way.  You can't remove them with a filter.

You're essentially sampling the pure mathematical signal without using an antialiasing filter first.  It's not possible to add the antialiasing filter after sampling.  You will just filter the already-aliased signal.  Bandlimiting has to occur before sampling.  You need to change the way you generate the waveforms if you want them to be bandlimited.  You could generate them using additive synthesis, as I've demonstrated here, or you could generate them the simple way at a much higher sampling rate to minimize the aliasing, then antialias filter, then decimate.  Or you could use BLIT or BLEP, etc.  (I think those are both "reduced-aliasing" not 100% eliminated, but I'm not sure.)


A question about spectral leakage I have read about spectral leakages in other posts here. From what i understand, it occurs because you dont have an integer number of time periods in your sampled data. By default a rectangular window (=no window) is used and it corresponds to a sinc function in the frequency domain. If we have an integer number of time periods in our data the zeros of the sinc cancel out every other frequency component giving us the correct value. If however we dont have an integer number of time periods in our sample, that corresponds to a shifted sinc and this is convolved with the frequency spectrum of the original resulting in leakage. In essence there is always leakage but it is masked by carefully choosing the number of samples to be an integer number of time-periods.

I wrote a simple program in scilab to see if i can remove this leakage with windowing. f(x) = Asin(20*pi*X). Maximum frequency = 10, sampling rate Fs = 30. I use a sinc window. A finite sinc window leaks in the frequency domain. For n samples the FFT(sinc) = (m1, m2, m3,...mn), m2 to mn are mirrored about the center. I changed the FFT(sinc) to FFTnew = (-m1, m2/abs(m2)-m2,-m3,-m4,...,mn/abs(mn)-mn) so that abs(FFT + FFTnew) = (0,1,0,...,1). I now create a new window called modifiedSincWindow = sinc + InvFFT(FFTnew). (in scilab invFFT(FFT(f(x)) = f(x), no scaling issues). the abs(FFT(modifiedSincWindow)) has all 0s except for 2 1s as expected. This is obvious from the linearity property of the DFT.

a) NumSamples = 30. FFT(f(x) x modifiedSincWindow) gives me the correct amplitude A.

b) NumSamples = 31. FFT(f(x) x modifiedSincWindow) does not give me the correct amplitude A. it doesnt even look anything like in (a). why is this? if we are convolving the frequency spectrum of the signal with the spectrum of the modifiedSincWindow and since the component falls in some bin shouldnt we get the same result as in (a). the spectrum for the modifiedSincWindow is exactly 1 for a particular bin and zero for the rest. And convolving this with the spectrum of f(x) should get me the same result as (a) but it isnt. Can anyone help me understand what exactly is happening here? why is it not working as expected? thanks for your comments.

EDIT: Added scilab code below.

My question basically boils down to if we are getting leakage because of the convolution by the sinc function, why not multiply the signal with a window whose magnitude response is 1 for one bin and zero for the rest. This should give us the right amplitude value irrespective of whether the number of samples is an integer number of time periods or not. right?

maxFreq = 10;
T = 3*maxFreq;

numSamples = T; %%(b) numSamples = T+1
timeRes = 1/T;

%%sampling the sinc function

sincVals = zeros(numSamples);

for i = 1:1:numSamples
    xval = (i - numSamples/2) * timeRes;
    if(xval == 0) %%then
        sincVals(i) = 1;
    else
        sincVals(i) = sinc(xval);
    end
end

fftVals = fft(sincVals);

modifiedFFTArray = zeros(numSamples);

for i = 1:1:numSamples
    modifiedFFTArray(i) = -fftVals(i);
end

modifiedFFTArray(2) = fftVals(2)/abs(fftVals(2)) - fftVals(2);
modifiedFFTArray(numSamples) = fftVals(numSamples)/abs(fftVals(numSamples)) - fftVals(numSamples);

invModified = ifft(modifiedFFTArray);

%%sincDiffs magnitude response is 0 for all bins except for 2 bins where it is 1

sincDiff = sincVals + (invModified);
signal = zeros(numSamples);

for i = 1:1:numSamples
    xval = i * timeRes;
    signal(i) = 1.5*cos(2*pi*maxFreq*xval) * (sincDiff(i));
end

binRes = T/numSamples;

%%disp(binRes)
binID = 1+floor(maxFreq/binRes); %%indexing is from 1 so add 1
%%disp(binID)
fftVals = fft(signal);

plot(abs(fftVals))

I see you've already accepted one of the answers, but you (or others) might be interested in: R. J. Webster, “Leakage Regulation in the Discrete Fourier Transform Spectrum,” Proc. IEEE, vol. 68, no. 10, Oct. 1980, pp. 1339-1341. -- Among the statements:

"A point insufficiently emphasized in the literature is that leakage in the DFT spectrum is a phenomenon of calculation.  The continuous distribution of leakage energy is determined only by the form of the window, and the DFT is a kind of lattice through which this continuous function is observed.  Because the lattice is fixed by convention, and the leakage distribution shifts according to the frequency of the signal, the discrete leakage components become dependent on the signal frequency."

(regarding the familiar 'picket-fence' effect): "The reader should note, importantly, that this effect is not a property of the harmonic nature of (the signal frequency fzero), as seems often suggested, e.g., [5], [6], but of the essentially arbitrary way G(f - (the signal frequency fzero)) is evaluated.  The apparently special relationship between harmonic frequencies and the rectangular window, in which the zeros of the window transform coincide precisely with the DFT sampling frequencies, is neither significant nor unique."

The author then goes on to show an example of using a shifted sum-of-cosines window centered on f = 10.25 that reduces the leakage to zero beyond bin 11 of his FFT when given an input signal of that frequency.


Is it accurate to say in digital audio that, when a fader is down, then its value is "-$\infty$"? Is it accurate to say in digital audio that, when a fader is down, then its value is "-$\infty$"?

What has confused me is, whether it's really $-\infty$ or some lower bound of the amplitude. it could be that the $-\infty$ is just a "joke" by the developers and it doesn't have mathematical validity.

If by "its" you mean the value of the fader, then yes - it's absolutely correct.
Fader defines the attenuation of the signal with respect to the reference level. The units are in logarithmic scale.

For example, no amplification (or in other words, multiplication of the signal by 1) corresponds to $0\; \mathrm{dB}$ since $10\log_{10} 1 = 0$.

If you want to completely attenuate your signal, then you have to multiply everything by $0$, hence $10\log_{10} 0 = -\infty$.


Can a Fourier Transform exist even if the j$\omega$ axis is not in the Region of Convergence in it's Laplace Transform A couple of confusions have been occurred. The Signal I'm considering is f(t) = sin(t)*u(t)

Fourier Transform of it can be derived.
$-i \pi  (\delta (\omega -1)-\delta (\omega +1))$
According to my mathematica code, the ROC of LaplaceTransformation didn't have the j$\omega$ axis in it's region of convergence.( Re{s}>0 )

So it's not stable. (Makes sense. sin(t) is not absolutely summable)
https://www.wolframcloud.com/objects/ramithuh/Published/misc_sin_laplace.nb

So, can a Fourier Transform exist even if the j$\omega$  axis is not in the Region of Convergence in it's Laplace Transform?

Things got worse when I considered f(t) = sin(t). It's Laplace Transform Integral didn't converge. So considering sin(t)*u(t) and sin(t)*u(-t) separately, I got two different ROC which doesn't have an overlap. Re{s} > 0 and Re{s} < 0. So it means that the Laplace transform of Sin(t) doesn't exist right? Initially what i thought is Laplace Transform can cover all the signals which FourierTransform covers.
Turns out it's not the case?

Please point out in which step of my reasoning is wrong...

Thanks a bunch! :)

Update: Thank you for pointing out my mistake. :D
$ \mathcal{F(sin(t)*u(t))} = -\frac{1}{2} i \pi  \delta (\omega -1)+\frac{1}{2} i \pi  \delta (\omega +1)-\frac{1}{2 (\omega -1)}+\frac{1}{2 (\omega +1)}$

You're right that the Laplace transform is not more general than the Fourier transform. They are just different. There are several (theoretically) important functions for which the Laplace transform doesn't exist, but the Fourier transform does. A few examples are

$x(t)=e^{j\omega_0t}$
$x(t)=\sin(\omega_0t+\phi)$
$x(t)= \textrm{sinc}(\omega_0t)$
$x(t)=\frac{1}{\pi t}$
$x(t)=\textrm{sign}(t)$

The first two involve Dirac impulses in their Fourier transform, whereas the third and the fourth have discontinuous Fourier transforms. I hope it is obvious why the complex exponential and the sine are important. The $\textrm{sinc}$ function is needed to represent the impulse response of ideal frequency-selective filters, such as low pass and high pass filters, and $1/(\pi t)$ is the impulse response of an ideal Hilbert transformer. Note that those functions are non-causal (i.e., they are not zero for $t<0$). The class of non-causal functions for which the (two-sided) Laplace transform exists is quite restricted because damping through multiplication with $e^{st}$ only works either for $t>0$ (if $\textrm{Re}\{s\}<0$) or for $t<0$ (if $\textrm{Re}\{s\}>0$), but not for both.

In your first example ($\sin(t)u(t)$) you got the Fourier transform wrong, even though it does contain Dirac impulses. The Fourier transform has Dirac impulses whenever the Laplace transform (if it exists) has poles on the imaginary axis.


calculating DC value from positive frequency with its DC undefined Is there a way of calculating DC value from just the positive frequency that has its DC undefined ? Simply if we remove the DC and the negative frequency components of the signal is there a way to reconstruct it from just the positive frequency values.

Assume you have the signal

$$x(t)=a+b\cos(\omega_0t)\tag{1}$$

with some non-zero real-valued constants $a$ and $b$. Now remove DC and the negative frequencies to obtain

$$x'(t)=\frac{b}{2}e^{j\omega_0t}\tag{2}$$

Can you calculate the value $a$ from $(2)$?


Proof of complex conjugate symmetry property of DFT According to the Proof :

  \begin{align}
X_n &= \sum_{k=0}^{N-1}x_ke^{-j\frac{2\pi k n}{N}}\\
X_{N-n} &= \sum_{k=0}^{N-1}x_ke^{-j\frac{2\pi k (N-n)}{N}}\\
&=\sum_{k=0}^{N-1}x_k e^{-j 2\pi k}e^{j\frac{2\pi k n}{N}}
\end{align}

  Using $\exp(-j2\pi k) = 1 \quad \forall \ k$

  \begin{align}
X_{N-n} &= \sum_{k=0}^{N-1}x_k e^{j\frac{2\pi k n}{N}}\\
\end{align}

How is $\exp(-j2\pi k) = 1 \quad \forall \ k$ true ? Doesn't it mean $-j2\pi k = 0$ ? But that's not possible right?

Hint:

According to Euler's formula we have
$$e^{-j2\pi k}=\cos(2\pi k)-j\sin(2\pi k)=\ldots$$


Phase Contribution of Zeros in Z-plane I have been studying Proakis, and stumbled upon minimum/maximum/mixed phase concept.

He said that all zeros inside unit circle in Z-plane contributes zero phase change, and all zeros outside unit circle contribute $\pi$ phase change. That is if $\Omega(\omega)$ is phase of the Z-transform, then, $\Omega(\pi) - \Omega(0)=\pi$.

Geometrically this does not make sense to me. Take the following example (The circle in both the pictures is an unit circle). This example is based upon the following transform :

$$H(z) = z^{-1}(z+ 1/2)$$

In the above picture, it is a single Zero system. The Zero is denoted by $Z(=-1/2)$ above. Here the $\Omega$ is the phase angle of the transform. Now say we rotate by more than $\pi/2$. And we get to the following picture,

Here $L_2$ is at say $\omega=3\pi/2$ angle. Then the new phase angle as $\omega $ approaches $\pi$ is going to be $\pi$. So from first picture at $\omega=0, \Omega=0$, and $\omega=\pi, \Omega=\pi$. So how come the phase difference $\Omega(\pi)-\Omega(0) = 0$ as claimed by the book for this transform. It also said that all zeros inside the unit circle would lead to zero phase change. But as described above, it does not make sense to me geometrically.

Can someone please guide me on where I am going wrong!

You're right that the angle $\Omega(\omega)$ in your sketch will equal $\pi$ for $\omega=\pi$. However, $\Omega(\omega)$ is not the phase angle of the system. In order to see this, let's write down the frequency response of a system with a single zero at $z_0=re^{j\theta}$:

$$H(e^{j\omega})=1-re^{j\theta}e^{-j\omega}=\frac{e^{j\omega}-re^{j\theta}}{e^{j\omega}}\tag{1}$$

From $(1)$, the phase angle is given by

$$\angle\{H(e^{j\omega})\}=\angle\{e^{j\omega}-re^{j\theta}\}-\angle\{e^{j\omega}\}=\Omega-\omega\tag{2}$$

So you indeed need the angle $\Omega$ as defined in your figures, but you also need to subtract the angular frequency $\omega$.

Let's assume $\theta=\pi$, as in your example. For $\omega=\pi$ (and $r<1$) you get $\Omega=\pi$, and, consequently, a phase of $\Omega-\omega=0$. Note that if $r>1$, you'll get $\Omega=0$ at $\omega=\pi$, and, consequently, a phase angle of $(-)\pi$.


Gaussian Pyramid - How is Subsampling Rate Related to Sigma? I found a gaussian pyramid implementation in a MOPS paper (feature detection). They use sampling rate $s=2$ and $\sigma=1$ - i.e. to generate a new level of the pyramid, the current level is smoothed with Gaussian blur of that sigma and then subsampled. The same parameters are used to build each new level.

However, I need to use smaller sampling rates, e.g. $s=1.5$, to get more pyramid levels (non-integer sampling rates would be achieved by interpolation).

How to choose appropriate sigma, knowing the sampling rate?

My first guess is to use $\sigma=\sqrt{s/2}$, since the variance of the gaussian filter is half the sampling rate (radius) and sigma (standard deviation) is square root of that quantity. But I am not sure if that's correct.

In another words: Given a sampling rate, I need to pick gaussian blur sigma preventing aliasing.

I've spent a bit of time searching for an answer to this as well.

I suspect since the Gaussian kernel has infinite support there will always be a little bit of aliasing, as there will always be frequencies passed that are greater than 0.25 of the sampling rate. That is, in the upper half of the spectrum.

So one would choose $\sigma$ to keep the amount of energy in the aliased bands to an acceptable level. Using $\sigma = 0.25 f_s$ would give 68%, $2\sigma = 0.25 f_s$ would give 95%, and so on.

From wikipedia:


How is Linear Canonical Transform a generalization of Fractional Fourier Transform? I have studied that Fourier transform changes the domain of a signal from time to frequency, and in that way it is a 90 degree shift. When it comes to Fractional Fourier Transform a generalization of Fourier Transform, the resulting transform can lie anywhere between time and frequency domain depending on parameter 'alpha'.

My question is how in the same way Linear Canonical Transform is a generalization of Fractional Fourier Transform.Its physical interpretation in terms of signal processing.
Can any body help?

The linear canonical transforms are all area preserving, orientation preserving, linear transforms on the time-frequency plane. The fractional Fourier transforms are a subset, namely the rotations. Both sets form (Lie-)groups under composition, and the rotations are a subgroup.


In case of Complex DFT spectrum, why the x axis range from N/2 to N point mean a negative frequency? When we transform  a complex signal into frequency signal by using a complex DFT,

The range from N/2 to N point on the X axis of the spectrum mean a negative frequency...

But i cant understand why it is..

In addition, we should not use the range over the N/2 point because of the aliasing issue.  But i heard that the negative frequency region is used in case of complex dft spectrum of a complex signal.

Very confused..

Could you clearly explain it for me?

Thank you.

it's because the DFT is periodic with period $N$.

$$  X[k+N] = X[k] \qquad \forall k \in \mathbb{Z} $$

this periodic extension is also true about the discrete-time signal $x[n]$

$$  x[n+N] = x[n] \qquad \forall n \in \mathbb{Z} $$

doesn't matter if either time-domain or frequency-domain representations are complex or not.

that means for $-\frac{N}{2}\le k < 0 $ (the negative frequencies), then $X[k] = X[k+N]$.  the explicit negative frequency component $X[k]$ has a copy $X[k+N]$, but that is in latter half of the DFT result.

If  $-\frac{N}{2}\le k < 0 $, then  $\frac{N}{2}\le k+N < N $.


How to calculate the luminance of an image? Is there any method for calculating luminance of an image?

As from wiki got to know Luminance is an amount of energy perceived by the observer and LUX is different from brightness.

But still is there any way to calculate LUX value of an image using any of image processing technique.

https://stackoverflow.com/a/4876528/4523099

  Just convert your image to YUV format and calculate the average of luma channel. Color conversion is a typical operation and any decent image processing framework supports it. For example, OpenCV has CvtColor.

However, to get an actual Lux value, you will need to calibrate your camera and tell whether it has an auto-setting aperture (see the link above for more details).


Remove unwanted object of image in MATLAB? I am trying to remove the unwanted object. I have an image that includes a house and the trees. I am trying to remove the trees. I am using 'imsubtract' to remove small edge. But the result is not good. The tree on wall still exist. How can I remove the tree in front of the building?

This is my code:

clear;

I = imread('Frame02.png');

BW = im2bw(rgb2gray(I), 0.5);figure;imshow(BW);

BW1 = BW;

CC = bwconncomp(BW);

numPixels = cellfun(@numel,CC.PixelIdxList);

[biggest,idx] = max(numPixels);

BW(CC.PixelIdxList{idx}) = 0;

figure, imshow(BW);

figure, imshow(BW1);

Ir = imsubtract(BW1,BW);

figure;imshow(Ir)

What you are asking is probably hopeless.

There are two quite difficult aspects:

detecting the trees that you want to remove. Given the very poor contrast, this will be hard to achieve by automated means. You need to implement some image segmentation method, then select the regions corresponding to the trees.
reconstructing the background image under the occluded parts. This is called image inpainting and requires sophisticated approaches.

All of this is highly technical and the standard tools in Matlab will probably not suffice.

If you just have one image to process, photoshopping will be much more efficient.


Combining 2 sinusoids of equal amplitude with different frequencies into 1 new "wave" This question is about combining 2 sinusoids with frequencies $\omega_1$ and $\omega_2$ into 1 "wave shape", where the frequency linearly changes from $\omega_1$ to $\omega_2$, and where the wave starts at phase = 0 radians (point A in the image), and ends back at the completion of the at $2\pi$ radians (point E), resulting in a shape similar to this, assuming $\omega_1$ is a lot smaller (larger period) than $\omega_2$:

Note that the amplitude is the same for both waves: points B and D are at equal distance from the $x$-axis.

Is there a closed-form definition for this, like linear frequency modulation?
Part of my question is also: where does point C end up on the x-axis if the "wave" starts at the origin (0,0)?

I was thinking about something like this:

1*sin((w1*t)+((w2-w1)/T)*((t^2)/2))

with:

w1: frequency 1
w2: frequency 2
t: a vector going from 0 to ((1/w1)+(1/w2))
T: some arbitrary number: this is part of the question I suppose..

Let us assume the following signal form $$x(t)=\sin(2\pi t f(t))$$ where $f(t)$ is some "instantaneuous frequency", which should increase linearly from $w_1$ to $w_2$ over the time.

The time should end, when the sine wave has done one period, i.e. its argument became $2\pi$. At this time, $f(t)=w_2$. From that we can calculate:

$$2\pi w_2 t_e=2\pi\\t_e=1/w_2,$$ where $t_e$ is the end time of the simulation. Hence, we have $$f(w)=w_1+\frac{w_2-w_1}{t_e}t=w_1+(w_2-w_1)w_2t.$$ We can plot this function:

w1 = 5
w2 = 10
T = 1/w2

Fs = 100000

t = np.arange(0, T, 1/Fs)

f = np.linspace(w1, w2, len(t))
plt.plot(t, np.sin(2*np.pi*f*t))
plt.title('Graph')
plt.grid(True)

Now, regarding the question, where the zero is. The sine goes to zero, when its argument becomes $\pi$. Hence, we have

$$2\pi (w_1+(w_2-w_1)w_2t)t=\pi,$$

which we can solve by (its a quadratic equation, so nothing special, but I use mathematica for convenience):

$$\text{Solve}\left[t (t \text{w2} (\text{w2}-\text{w1})+\text{w1})=\frac{1}{2},t\right]$$

$$\left\{\left\{t\to \frac{\sqrt{\text{w1}^2-2 \text{w1} \text{w2}+2 \text{w2}^2}-\text{w1}}{2 \left(\text{w2}^2-\text{w1} \text{w2}\right)}\right\},\left\{t\to \frac{\sqrt{\text{w1}^2-2 \text{w1} \text{w2}+2 \text{w2}^2}+\text{w1}}{2 \left(\text{w1} \text{w2}-\text{w2}^2\right)}\right\}\right\}$$

There are two solutions, because for negative t, we can also end up with a solution (both $f(t)$ and $t$ would be negative to yield 1/2 as the argument).

To answer the comment from @Jason R: This indeed can also be seen as a frequency-modulation, given by

$$x(t)=\sin\left(2\pi\int_0^t\phi(\tau)d\tau\right)$$

when setting

$$\phi(\tau)=w_1+2(w_2-w_1)w_2\tau.$$

In this case we have $tf(t)=\int_0^t\phi(\tau)d\tau$.


Random Noise Removal In Images I'm trying to remove random noise from an image. The noise is arbitrarily random so it must not necessarily follow a certain distribution (Not every pixel has the same chance of being affected by noise and not every pixel is affected by the same amount of noise). Such an image might for example look like this example image:

The usual fast'n'cheap filters (like average etc.) just do not cut it for such images. It's pretty obvious that there's still enough information in the image left for it to be reconstructible but I don't think there's yet a tool available that can actually do that. Any hints on how I could try reconstruction? Tools, algorithms?

Update 1

The closest I've come so far is this:

It's far from being good but just to show what I mean by removing the noise.
The resulting image should have more or less uniform areas of about the correct colors.

An ideal outcome would roughly look like this:

Update 2

Example with little noise

Example with more noise

Update 3

PNG images:

For the lena fans out there:

There are a lot of ways to do image denoising.

I will give several known examples:

Hard/soft thrsholding of the wavelet decomposition: Do a wavelet decomposition of the image and for zero all coefficients less than some minimum value and linear the other values.

Non Local Means: We use the fact that for "natural images" for each patch in the image we can find very similar other patches in the image. So for each pixel we look at the neighborhood around him and we search for other pixels that has the same neighborhood and than we use those patches to calculate the value of the original pixel.using the fact that the patches are very coorelated but the noise isn't.

Block Match 3D: We do the same as before finding similar patches around the image and than we arrange them as a 3D block and apply wavelet thresholding. This is more complicated and you can read about this. A basic scheme for this is:

There a lot more ways and just by googling it you will get a lot of information hope this basic intro will help you enough.


Does the system function H(z) for a filter always have symmetry above and below the real x-axis? When I look at the Pole/Zero diagrams in the DSP books it appears that whatever
pole or zero is above the real x-axis from:

$$\omega=0 \rightarrow \pi$$

is always reflected as a mirror image below the x-axis from:

$$\omega=0 \rightarrow -\pi$$.

(And whatever pole or zero is on the real axis with $Im\{H(z)\}=0$ is not reflected anywhere.)

For instance, if you have a zero at $e^{j\omega_0}$ and ($\omega_0\ \ne \pi$, and  $\omega_0\ \ne 0$) then you are always going to get a zero that is below the x-axis at $e^{-j\omega_0}$ regardless what type of filter or function the system function H(z) is implementing.

I'm wondering what law or rule this relationship comes from?  and are there any cases where this is not true for a realizable system function of H(z)?

No it does not (always posses symmetry).

However, for real systems whose coefficients $a_k ~,~ b_k$ of system function,

$$H(z) = \sum_{n} h[n] z^{-n} = \frac{\sum_{k=0}^{M}b_k z^{-k}}{\sum_{k=0}^{N} a_k z^{-k}} $$ are real, then accordingly its roots (which are the poles and zeros of $H(z)$) are either real or complex-conjugates. And because of this conjugation, the poles (or zeros) on the Z-plane posses the vertical symmetry as you've recognised.

On the contrary, for complex systems whose system function $H(z)$ is a complex polynomial of complex coefficients, we cannot constraint the roots of denominator and numerator of $H(z)$ to be real or complex-conjugate, and hence they can posses non-symmetric poles and zeros...


How to decide the number of cepstral coefficients? For audio classification we are using cepstral coefficients, but how should I decide the number of coefficients?

From CMUSphinx FAQ:

  There are various types of MFCC which differ by number of parameters,
  but not really different for accuracy (it might be a few percent worse
  or better).

  The interpretation of MFCC (Roughtly introduced Alan V. Oppenheim and
  Ronald W. Schafer. From Frequency to Quefrency: A History of the
  Cepstrum. IEEE SIGNAL PROCESSING MAGAZINE) is not applicable as such,
  and the use of 12 or 13 coefficients seem to be due to historical
  reasons in many of the reported cases. The choice of the number of
  MFCCs to include in an ASR system is largely empirical. To understand
  why any specific number of cepstral coefficients is used, you could do
  worse than look at very early (pre-HMM) papers. When using DTW using
  Euclidean or even Mahalanobis distances, it quickly became apparent
  that the very high cepstral coefficients were not helpful for
  recognition, and to a lesser extent, neither were the very low ones.
  The most common solution was to “lifter” the MFCCs - i.e. apply a
  weighting function to them to emphasise the mid-range coefficients.
  These liftering functions were “optimised” by a number of researchers,
  but they almost always ended up being close to zero by the time you
  got to the 12th coefficient.

  In practice, the optimal number of coefficients depends on the
  quantity of training data, the details of the training algorithm (in
  particular how well the PDFs can be modelled as the dimensionality of
  the feature space increases), the number of Gaussian mixtures in the
  HMMs, the speaker and background noise characteristics, and sometimes
  the available computing resources.


Why is scale space (DoG) needed to detect scale invariant features? I would have a theoretical question, about the detector part of the SIFT algorithm, which as it is explained by D. Lowe in his paper, used the DoG to detect keypoints in an input image.

Problem formulation: Let $x$ be my original image, and let $y$ be a scaled-by-two version of $x$.

In the SIFT algorithm: I don't understand how the SIFT detector can detect a same given keypoint from image $x$ in image $y$. I don't understand how the "scale space" DoG pyramid is used exactly to achieve this "scale invariance".

Please note that I have watched quite a few videos and read a few tutorials about it. I understand how the DoG pyramid is constructed, by successive blurring of the image. I also understand that keypoints are the maxima or minima from comparing to adjacent images in the DoG pyramid. (I understand that the difference of 2 blurred images with different variances gives an edge-like image.)

However, for me, in every tutorial, they sort of "skip" the fundamental part that would clearly explain why it makes it scale invariant.

In other words, I understand that different gaussian blurs (increasing variance) allow to highlight different types of details in images (blob detector). I'm ok with that, but I don't see any relation or reason intuitively why this has anything to do whith "scale" and why and how that is used to detect a same keypoint in image $x$ and $y$ with for example image $y$ zoomed in twice.

I have essentially the same (unanswered) question as him.

As you know, you find interest point (SIFT points) by finding the local maxima in scale-space, it mean the response of the detector must be maximum regarding the coordinates and also the scale. So knowing the scale you select features from that scale. Now consider the same process for same interest point in the same picture but with different zoom, the detector maxima will be in the same scale, so it will collect the same feature. So in the matching these two same feature will be matched together.

So the scale invariant is a result of finding maxima with respect to scale.


How to solve equations for the MMSE receiver? Given the following equations for the achievable rate of the Minimum Mean-Squared Error Receiver [1]:

$$\mu = \frac{1}{K-1} \sum_{i=1, i \neq k}^{K}{\frac{1}{Mpd_{i}\left(1 - \frac{K-1}{M}+ \frac{K-1}{M}\mu \right) + 1}}$$

and

$$\kappa \left(1 + \sum_{i=1, i \neq k}^{K}{\frac{pd_{i}}{ \left(  Mpd_{i}\left(1 - \frac{K-1}{M}+ \frac{K-1}{M}\mu \right) + 1 \right)^2} } \right) = \sum_{i=1, i \neq k}^{K}{\frac{pd_{i}\mu + 1}{ \left(  Mpd_{i}\left(1 - \frac{K-1}{M}+ \frac{K-1}{M}\mu \right) + 1 \right)^2} }$$

I've got to find, $\mu$ and $\kappa$. I'm trying to solve these equations with matlab, firstly I tried with solve and then vpasolve, however, when $d$ is a vector with different numbers (i.e., random) it seems matlab never finds a solution. For the case when $d$ is a vector with all values equal to 1 for example, it finds the solution quickly.

This is the script I'm using to find the solution:
https://pastebin.com/AWDwyR0U

Am I doing something wrong? Is there a faster way to find a solution to these equations?

Reference:

[1] Hien Quoc Ngo, Erik G. Larsson, and Thomas L. Marzetta, "Energy and Spectral Efficiency of Very Large Multiuser MIMO Systems", IEEE Transactions on Communications, vol. 61, no. 4, April 2013.

I doubt there is general closed-form solution. For the first equation, you essentially want to solve $$\mu=\sum_k \frac{1}{a_k+b_k \mu}$$ for $\mu$ which equals to search for the zeros of a $k+1$th order polynomial, which is impossible for $k>3$. So, I fear you need to do a numeric solution. Once you have $\mu$ you can go for a simple calculation of $\kappa$.

Regarding a numerical solution, here's a plot of $$\mu-\sum_k \frac{1}{a_k+b_k \mu}$$ which should become $0$ for the solution:

N = 7
d = np.random.rand(N)
K = N+1
M = 16
p = 1

f = lambda mu: mu - 1/(K-1)*sum(1/(M*p*d_i*(1-(K-1)/M+(K-1)/M*mu)+1) for d_i in d)

mu = np.linspace(-5,5,10000)
plt.grid(True)

plt.plot(mu, f(mu))
plt.ylim((-0.5,0.5))

So, there are multiple solutions. Maybe you need to pick to positive one?


Spectral location of individual OFDM symbols The generic IFFT OFDM transmitter takes $N$ serial input data symbols, parallelizes them, maps them to their signal constellation points, and processes all $N$ resultant complex values using the IFFT to produce $N$ subcarriers (assuming no virtual carriers).

The IFFT takes $N$ input values and produces $N$ output values, such that changing a single input value will change all output values. Therefore each input symbol contributes to each of the $N$ subcarriers; each symbol is "shredded" across all subcarriers.

However, several sources state that in OFDM "each modulation symbol is mainly confined to a relatively narrow bandwidth" (Communications Engineering Desk Reference, Elsevier, 2009, 1st ed, p.287).

Are these two different OFDM implementations, or is the IFFT implemented in a different way to achieve this?

The IFFT transforms a discrete signal to the time domain. Therefore, the input signal of the IFFT should be interpreted as a frequency domain signal.

For OFDM this means that the subcarriers are defined before the IFFT. The IFFT ouput is the time domain representation that will eventually be transmitted over the channel. Before it can be transmitted it must first be converted to an analog signal. Analyzing this signal in time domain, e.g. by viewing it on an oscilloscope, won't tell you anything about the subcarrier bandwidth.

To see the bandwidth of one OFDM subcarrier you have to analyze the analog signal in the frequency domain, e.g. by using a spectrum analyzer. If you modulate all subcarriers the spectrum will look nearly rectangular, with a flat top, because the subcarriers are close together and overlap. If you modulate only one subcarrier (all inputs of the IFFT are set to zero except one), the spectrum will be much narrower, its bandwidth approximately reduced by a factor that is equal to the number of subcarriers. And it will have a sinc-like shape.


Image processing: using PCA for 3d object recognition I read this lecture on using PCA/eigenfaces: http://www.cse.psu.edu/~rcollins/CSE486/lecture32.pdf

And I want to use this general procedure to recognize a 3D clay model that I took pictures of from multiple angles.

I'm having trouble conceptualizing how it's actually implemented in software. I've "trained" the initial mean and eigenvectors for using ~40 images offline. So I have new input images that will have my 3D model in it (or wont).

My confusion: I'm worried that there's so much else going on in my input image that applying the projection + thresholding equation isn't going to work. It almost seems like I need to sub-image my input many times and try over and over again to get the right projection (almost like multi-template matching). Or maybe I can detect blobs first, create sub-images from the blobs, then try the projection.

Because I have multiple angles of the model I suppose that handles rotational invariance, but what about scale invariance? Seems like whatever subimage I get from my input, I have to try at different scales. Or I'm not understanding something fundamental about how to use PCA.

Thank you for any help; for reference I would like to implement this in OpenCV (Java-based).

Principle Component Analysis by itself is not going to provide rotational or scale invariance. You still have to provide templates for rotational and scaled models for image/template matching. What PCA does do is reduce the dimensionality of the search space  by converting a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables called principal components.

Now since number of principle components <= number of components, you save on the search space.

What the lecture tries to tell you is with the eigenvectors, you get an efficient representation of the scales of the template which are easier to manipulate(a single scale factor) then manipulating the entire image everytime.


Discrete Fourier Transform - measuring signal at equi-spaced points I found a really intuitive explanation of Discrete Fourier Transform here.

To summarize the author's approach briefly :

What if any signal could be made from circular motion (repeating cycles) ? This is the central theme of Fourier Transform
Every circular path needs a size, speed, and starting angle (amplitude/frequency/phase).
The combined "position" of all the cycles is our signal.

The author provides an interactive visualizer (screenshot below):
The cycles are shown on the left and the paths shown on the right. It shows 3 cycles (a cycle of 0 Hz, a cycle of 1 Hz and a cycle of 2 Hz). The paths of the cycles are in green. The combined path (signal) is in blue.

The author says

  The yellow dots are when we actually measure the signal. With 3 cycles defined (0Hz, 1Hz, 2Hz), each dot is 1/3 of the way through the signal. In this case, cycles [0 1 1] generate the time values [2 -1 -1], which starts at the max (2) and dips low (-1).

Why are the dots spaced so ? More generally, if there are $n$ cycles, why is the signal sampled at $n$ points equi-spaced around the circle ?

The "Cycles" entries are the amplitudes of cosine waves of frequency 0Hz, 1Hz and 2Hz respectively, so:

0Hz:  $0\cdot\cos(0)$,

1Hz:  $1\cdot\cos(2 \pi t)$,

2Hz:  $1\cdot\cos(4 \pi t)$

The time values are amplitudes of the resulting sum at specific sample points in time.
The time runs over 1 cycle of what in Fourier series analysis is call the fundamental frequency which in this case is 1Hz.

The Discrete Fourier transform is really a Fourier series, so the signal being represented is always periodic and always has a maximum time equal to 1 period of the fundamental frequency which represents one trip around the circle.  Every other frequency represented in the Fourier transform is an integer multiple of the Fundamental frequency or represents an integer multiple of trips around the circle.

When computing a discrete Fourier transform a signal is first sampled at some rate that must be higher than the frequencies represented by the transform.  For traditional frequency analysis, sampling is always done using uniform spacing.  This is called uniform sampling and since it must happen faster than the signal you are representing, you must have more than 1 sample uniformly spaced around your circle.  To be more precise, you must sample your signal at least 2 times higher than the highest frequency that you want to represent in your analysis, so for this case the sampling should be a minimum of 4Hz.

A sampling rate of 1 Hz would give you a single yellow dot on the circle at 0 angle (x=1).
2 Hz gives you 2 yellow dots, one at 0 angle and one at $\pi$ (x=1 and x=-1).
In this example you should have a minimum of 4 samples to represent the frequencies illustrated, so you should have 4 yellow dots equally spaced around the circle starting at 0.

The circle is usually represented on a complex plain with the vertical axis representing imaginary numbers.  If this is done, the corresponding dots correspond to specific times calculated from the complex exponential $\exp(i 2 \pi t)$.  With the sample times starting at t=0 and progressing as 0,T/4,T/2,3T/4 (T being the period of the fundamental), the dots correspond to placement around the circle at the angles $0, \frac{\pi}{4}, \frac{\pi}{2},$ and $\frac{3 \pi}{4}$.


pre-emphasis filter - parameter a Can someone give my a link to study to solve this problem.
I can not find any formula for this problem...

Suppose we filter singals with the below pre-emphasis filter:

$y(n)=x(n) - ax(n-1)$

How can I compute $a$ if I know that:

fs=10.000Hz
length 100 in which 98 is zero.
Frequency response at 2100Hz is 1.3429 dB

I don't want the solution only a hint I found this in my notes: 10log10(a^2*w^2+1)=db but with this formula i don't use some information which given.

Thanks

Hope this is right...
$$h[n]=\delta[n]-a\delta[n-1]$$
$$H(e^{jw})= 1- ae^{-jw} =(1-a\cos(\omega))+ja\sin(\omega)$$
$$|H(e^{j\omega})|= \sqrt{1 + a^2 -2a\cos(\omega)}$$
$$DB_\omega \triangleq 20\log(|H(e^{j\omega})|)$$
$$20\log((1 + a^2 -2a\cos(\omega))^{1/2})= DB_\omega$$
$$1 + a^2 -2a\cos(\omega)= 10^{DB_\omega/10}$$
$$a^2 -2a\cos(\omega) +1-10^{DB_\omega/10}=0$$

with $\omega = 2 \pi \frac{2100}{10000}, \ DB_\omega=1.343 \ \Rightarrow \ a \in \{0.9, -0.4026\}$

$a = -0.4026$ results in a low-pass pre-emphasis whereas $a = 0.9$ results in a high-pass pre-emphasis.


Filter impulse response order calculation I want to generate the following periodic sequence:
h[n]={1,2,3,4,0,0,0,0,1,2,3,4,0,0,0,0,1,2...}.
I am planning to generate this as an impulse response of a filter with an order no more than 8.

I start by finding the z transform of h[n] but this is infinite sequence and can't achieve an order less than 8.

Any hint will be appreciated.

You need a system with order $8$, because the period of the desired sequence is $8$. Obviously, you'll need a recursive system to generate such an infinite sequence. Since this looks like another homework problem, I won't write out the complete solution but I'll give you a few hints:

Let $g[n]=\{1,2,3,4,0,0,0,0\}$ be the first period of $h[n]$. Now note that

$$g[n]=h[n]-h[n-8]\tag{1}$$

The $\mathcal{Z}$-transform of (1) is

$$G(z)=H(z)(1-z^{-8})\tag{2}$$

from which

$$H(z)=\frac{G(z)}{1-z^{-8}}\tag{3}$$

follows. Now it only remains to find the $\mathcal{Z}$-transform $G(z)$ of $g[n]$. Then you can use Eq. (3) to compute $H(z)$ from $G(z)$. If you realize a filter with transfer function $H(z)$ (which is a recursive filter of order $8$), then its response to a unit impulse will equal the desired sequence $h[n]$.


Studying spectral leakage due to DFT length and frequency resolution I guess this is a pretty basic question, but I am kinda stuck... So if I am remembering the DSP theory correctly, spectral leakage can occur when we take "inappropriate" combinations of the following:

1) FFT length

2) sampling frequency

3) fundamental frequency of input signal

Here is a MATLAB script for demonstrating this, in a case of a sinusoid input:

F0=10; %Frequency of the sinusoid
Fs=100; %Sampling Frequency
observationTime = 1; %observation time in seconds
t=0:1/Fs:observationTime-1/Fs; %time base

x=sin(2*pi*F0*t);%sampled sine wave

N1=100; %DFT length same as signal length
X1 = 1/N1*fftshift(fft(x,N1));%N-point complex DFT of x
f1=(-N1/2:1:N1/2-1)*Fs/N1; %frequencies on x-axis
stem(f1,abs(X1));

MY thought was that if the fundamental frequency is a multiple of the frequency resolution, then there will be no spectral leakage, since there is a bin that corresponds to the exact value of the fundamental frequency F0. However, it seems like I am wrong here...

Can someone please explain to me what's the flaw in my rationale? Why is there leakage in the 3rd plot? Thanks in advance

It looks like what you are doing is only changing the DFT size (using N1 variable), this is equivalent to padding your signals with zeros up to a given length. Your signal (for N1=200) looks like:

It can be immediately seen that the periodic extension is not really a continuous waveform. This discontinuity is introducing the leakage effect that you are observing.

Below a quick and dirty code that demonstrates that. I did both padding of input signal in time domain, as well as calling the fft function with twice the signal size. Plots are in dB scale.

F0=10; %Frequency of the sinusoid
Fs=100; %Sampling Frequency
observationTime = 1; %observation time in seconds
t=0:1/Fs:observationTime-1/Fs; %time base

%% Ideal case
x=sin(2*pi*F0*t);%sampled sine wave

N1=100; %DFT length same as signal length
X1 = abs(1/N1*fftshift(fft(x,N1)));%N-point complex DFT of x
f1=(-N1/2:1:N1/2-1)*Fs/N1; %frequencies on x-axis
X1 = 10*log10(X1/max(X1));
subplot(2,1,1)
plot(f1,X1, 'o-');
title('Ideal case')
grid on

%% Padding signal with zeros
x2 = [x zeros(1, length(x))];

N2=200; %DFT length same as signal length
X2 = abs(1/N2*fftshift(fft(x, N2)));%N-point complex DFT of x
f2=(-N2/2:1:N2/2-1)*Fs/N2; %frequencies on x-axis
X2 = 10*log10(X2/max(X2));
subplot(2,1,2)
plot(f2,X2, 'o-');
hold on
grid on

%% Padding with zeros implicitly
N2=200; %DFT length same as signal length
X3 = abs(1/N2*fftshift(fft(x, N2)));%N-point complex DFT of x
f3=(-N2/2:1:N2/2-1)*Fs/N2; %frequencies on x-axis
X3 = 10*log10(X3/max(X3));
plot(f3,X3, 'r-');
title('Padding with zeros')
legend({'Padded vector', 'Twice the FFT size'})


The impact on the frequency of adding zero samples and non zero samples Can someone give me a brief explanation of what is going on in here? I struggle to understand the full differences between the impact on the frequency of added zero samples and adding samples which are not zeroes. (In between the samples).
I learnt this course 100% from the lectures so if I have some misconception regarding the basics, I'd be glad if you could kindly help me understand.
The key points I know so far:

Adding zero(ed) samples in between the existing ones in the time domain, squizes the frequency domain and adds height to it.
The squized frequency domain is squized while it's normalised (once we turn it to be periodic with period of 2π).
Adding non zero(ed) samples to our time domain, we get an added zeroes vector to the end of the frequency domain's period. (Padding zeroes in frequency). Which means that when it's being normalized into 2π period's length, then it's also being squized.

From that we conclude, that adding zero samples and non zero samples do the same exact thing.
[Question:] where is my mistake in the above logic sequence? Where is the difference between the impact of adding zero samples in the time domain on the frequency vs adding non zero samples in the time domain on the frequency?

As soon as you sample in the nuiqist rate, sampling faster won’t add information(!). Sampling zeroes can be looked at as approximation mistake(!) It gives the same frequency as the original signal mathematically because it’s multiplying by zero and IF by contradiction, there would have been any other value instead of the original value, which WOULD [again by contradiction] add some information aside from what we have with the niquist samples, THEN we would need a higher niquist rate and that would be a contradiction to our given Max. Niquist frequency.
Yes, indeed sampled zeroes gives the same signal and they are considered as a measure mistake given the Niquist frequency!


selecting bricks from image I have an image of bricks of different sizes. What I want to do is when I click on the corner of the brick it should mark the lines around the brick. Can anyone please help me how can I do that...
here is the image

Let me suggest a simple solution , however if you want to go for a more robust solution let me know.

So first I would look for a way to enhance the signal. In my experience red objects appear even stronger if you look at red-green

diff = im(:,:,1) - im(:,:,2);
figure;imshow(diff,[]);

Than I would use simple thresholding (choose threshold value manually ) to get a blackwhite image.
To make it more robust I would use some morphology operations like opening and filling. (If you dont know what that is I would look at  matlab help.
Than finally you can select the brick with bwselect and take its perimeter with bwpeim.

Viola:

diff = im(:,:,1) - im(:,:,2) ;
bw = diff > 35;
bw = imopen( bw , ones(1,5) );
bw = imopen( bw , ones(5,1) );
bw = imfill(bw , 'holes');
brick = bwselect(bw);
brick = bwperim(brick);

and you can see a simple output here:
Please note that this is quite specific for this image but that's a good simple way to get you started.