What are the differences between the DWT and the MODWT? What are the differences between the DWT (Discrete Wavelet Transform), which is the most classical algorithm and the Maximum Overlap Discrete Wavelet Trasnform (MODWT)?

Both these algorithm are available in MATLAB. As far as I know MODTW is more suite for multi-resolution analysis, whereas DWT is more useful for signal de-noising and compression.

Thanks for your clarification.

from "wavelet methods for time series analysis"


Why there is static noise and voice distortion after convolving it with impulse response? I have been working on the decorrelation of audio signals while following this paper :
The Decorrelation of Audio Signals and Its Impact on Spatial Imagery
So far I have generated impulse responses in such a way that magnitude is 1 in all frequency spectrum, while random phase generation between -π to π. I correlated these signals and I am getting the required correlation.
Now moving on, when I convolve this impulse response with my audio signal, results in the case of decorrelation are good but the audio becomes noisy and distorted. I am attaching the reference wav file and decorrelated wav file.
Some points regarding computation:

All audio data computations are in float range [-1.0:1.0].
I am using Python to perform all the calculations
The input reference file is in Short (PCM 16 bit) format which is converted to float range [-1.0:1.0].

I have saved convolved_decorrelated_output.wav in both formats PCM 16 bit short and Float as mentioned above. Still getting the same results. So there are no issues regarding data type conversion as per my knowledge.
Any help will be appreciated.
Regards,
Khubaib Ahmad
P.S: Link to download wav files:
https://www.researchgate.net/post/Why_there_is_static_noise_and_voice_distortion_after_convolving_it_with_impulse_response
Impulse Response:

Convolution_code:
resultant_chunk = list()
for chunk in data_chunks:
    resultant_chunk.append(np.convolve(chunk,impulse1,mode="same"))

resultant_chunk = np.hstack(resultant_chunk)
# diff = len(resultant_chunk) - len(data)|
diff = len(resultant_chunk) - len(recorded_data)
resultant_chunk = resultant_chunk[:-diff]
len(resultant_chunk)

I came up with the solution after careful analysis of all the steps I was performing.
So basically, I was performing convolution on a 21ms chunk of data and the convolution results were different but when I convolved the complete audio signal with my impulse response, there was no distortion.
So, I had to implement the Overlap-add method to counter this issue.
Regards,
Khubaib


Fractional Resampler - GNU Radio Fractional Resampler block in GNU Radio has parameter phase shift but there is a lack of information in which factors user has to provide this parameter (radians, degrees or samples?). There is only an error message providing and info that the range of this parameter should be set between 0 and 1. I have created simple simulation with the fractional resampler block to kinda 'debug' its behavior. There is an 100 length input vector consisting of zeros but one in the middle (50th index). Then there is an interpolatior with the RRC filter. Given 0 phase shift and 1 for resampling ratio the output from the interpolating FIR Filter and Fractional Resampler are not in line, meaning that original signal is delayed by 3 samples and to compensate this phenomena I needed to delay samples by 3 in the branch where Fractional Resampler exists. Question are:

Why 3 samples offset occurs?
How to interpret phase shift parameter from Fractional Resampler block?

I don't know if its relevant information but QT Time Sink has normal trigger set to 0.1.
I am enclosing the flowgraph and the source
https://gist.github.com/marcinsztajn/afe03d05653864b10125fefe7d8b888d

"Phase shift" is poor terminology as time delay is not exactly the same, but the significance of the 0 to 1 factor is a fractional delay in units of samples.
The filter itself will always have additional delay, as a causal filter which is simply a processing delay independent of resampled delay as given by the phase shift parameter. For example a resampling delay of 0.5 will resample the original waveform such as to provide all new samples that are midway in between the current samples, and this new resampled waveform will be completely delayed by the filtering delay in time related to the input waveform.
The example in the figure below shows the original samples as blue dots of an underlying continuous time waveform as the line. If we used a "phase shift" parameter of 0.2, the samples as given by the red dots would be determined as the interpolated fractional delay values. These red dots would all be shifted even further to the right relative the original locations as shown in the lower plot due to the processing delay of the filter itself. (As we would see with any causal filter).


Use of pseudo-random binary sequence (Gold, Kasami, etc.) in frequency hopped code-division multiple access (FH/CDMA) I have a clear idea of how a PN-sequence is used in a Direct Sequence Code-Division Multiple Access(DS-CDMA) system with BPSK. In such a system, a different pn-sequence is assigned to each user and each transmitted bit is XORed with the binary pn-sequence.

Now my question is:

For each symbol time, how are this binary pn-sequences used in a FH/CDMA system to chose which hop frequency to use if there are more than two possible hop frequencies?

If I have 2 possible hop frequencies to use, the answer is trivial, because each bit of the pseudo-random sequence would tell which of the two hop frequencies must be used. But when having more than two frequencies... How is this selection done?

There are many ways to do this. A simple method is to take the current value of the register (not just its output) as the carrier frequency selector. For instance, a (6,1) LFSR with 63 states can select between 63 different channels.

Other methods are based on cubic and quadratic congruence codes. I don't know much about them myself, but look up the papers by Titlebaum et al.


How to change the angle unit labels on polar plot in MATLAB? I hope this is not off-topic. I'm wondering is there anyway to change the default angle unit label in MATLAB polar plot. The default angle label is from 0 to 360, at 30 increment. It goes like 0,30,60,... But I want to change it, say, from 0 to -360, at -30 increment. So the labels will be like 0,-30,-60,... Is there anyway to do it?

When you opens the polar plot in the figure editor you can find out that the coordinates around the plot are not from the axis.
It are hidden text objects (Press on one you see it). So this means you have to edit them by yourself. a work around for this is this:

figure
t = 0:.01:2*pi;
polar(t,sin(2*t).*cos(2*t),'--r')

hHiddenText = findall(gca,'type','text');
Angles = 0 : 30 : 330;
hObjToDelete = zeros( length(Angles)-4, 1 );
k = 0;
for ang = Angles
   hObj = findall(hHiddenText,'string',num2str(ang));
   switch ang
   case 0
      set(hObj,'string','0');
   case 30
      set(hObj,'string','-30');
   case 60
      set(hObj,'string','-60');
   case 90
      set(hObj,'string','-90');
   case 120
      set(hObj,'string','-120');
   case 150
      set(hObj,'string','-150');
   case 180
      set(hObj,'string','-180');
   case 210
      set(hObj,'string','-210');
   case 240
      set(hObj,'string','-240');
   case 270
      set(hObj,'string','-270');
   case 300
      set(hObj,'string','-300');
   case 330
      set(hObj,'string','-330');
   otherwise
      k = k + 1;
      hObjToDelete(k) = hObj;
   end
end
delete( hObjToDelete(hObjToDelete~=0) );

If you want other directions etc i think you do that manually otherwise the plot start tot mess up.


image reconstruction using a set of kernels I try to be brief. For what I understood Convolutional Neural Network CNN for style-transfer extract/learn the main features of source Image. I think these features basically are Kernels (of different sizes depending on Levels). What I ask is how to reconstruct image once we obtained the Kernels. To simplify let's assume I have N-kernels of the same size. How should I reconstruct the Image using only the Kernels ? If question is not clear I'll try to clarify it. (I try using and weighting their responses, but not succeed) Thanks.

PS I wrongly posted to stackoverflow the question... but I think this DSP is the right place. link to Stackoverflow: https://stackoverflow.com/questions/41833641/image-reconstruction-using-kernels

I think what you would like to do is something similar to an Auto-Encoder, in which one would append an Encoder to reconstruct back the image. The 'kernel', or more appropriately the intermediary representation is used to synthesize an image, and an appropriate reconstruction objective is minimized to recover the input. On the other hand, you could always train a second upsampling network, which reconstructs the input image, given a 1D/2D feature.


Why is my resonance phase lag discontinous? I am programming a phase-locked loop to track the resonance of a tuning fork so that the tuning fork is always on resonance. This will be used in condensed matter physics through atomic force microscopy.
Now, I am not certain about how to carry this out, so I was wondering if someone here could help me.
I have programmed an FPGA to output an arbitrary driving voltage, and this drives a tuning fork. This tuning fork then generates a signal that is eventually converted to a voltage and read by my device. Now, I also have a lock-in amplifier to detect the phase lag between my DAC (technically, the NCO that outputs via DAC) and ADC signals. I have made a graph of the phase difference after taking 4000 steps around the resonance frequency. Also, there is a graph of the LIA's detected magnitude*2.
I do not know what any of the transfer functions are, I apologize.

Phase difference (phase lag) graph on left. Resonance centered at 32.7635 kHz

However, my phase difference graph is likely incorrect. I found the following that displays continuous "pretty" graph that may be what I am supposed to have.
Phase lag in blue
Does anyone know why phase graph has a max/min of 1/-1? If so, does anyone know how I may fix it?
Thank you

You seem to have 3 problems
1 - You have a phase wrap. It happens when the phase goes past 180 degrees. The phase will wrap back to -180 degrees. You can fix this by unwrapping the phase.
2 - Your phase is in normalized radians instead of being in radians or degrees. This is not a problem per se but when you try to unwrap the phase, you should be aware that phase unwrap blocks usually expect a phase in radians, not in normalized radians.
3 - The measured phase seems to be inverted compared to what you expect.


System Identification using LMS Adaptive Filter I just have a question about using an least-mean-squares algorithim adaptive filter for system identification. Consider the following

I am told that as the error converges to a small value, the adaptive filter coefficients w[k] will indeed repersent the unknown system h[k]. Now, that doesn't make sense to me since the n[k] is being used in the error calculation.

Won't the adaptive filter coefficients w[k] now repersent the unknown system coefficients h[k] AND the noise n[k]?

Because an LMS estimator will, over time, "average out" uncorrelated zero-mean noise. It's pretty much in the name.

But yes, you're right, there is a noise component in the estimate; one of the qualities of an estimator is how little the noise variance influences the estimate variance after a given length of observation.

That's the case, however, for all estimators: you measure signal + noise, and you estimate parameters from that. The parameters must be somewhat noisy; otherwise, there's something broken with your noise model.


QPSK Phase Recovery : Maximum-likelihood Vs Costa's Loop I have implemented a QPSK Costa's loop in GNU Radio and run some BER benchmarking. All is going well. Out of curiosity, I'm trying to implement the Maximum Likelihood phase recovery shown below.

For the Costa's loop, a loop bandwidth, which specifies gain constants for the loop filter (F(z)) and the integrator,  of a few percents of the symbol rate is adequate.
My question is, first of all, since the ML phase recovery doesn't seem to have any traditional loop filter (e.g. proportional-plus-integrator), does the very concept of loop bandwidth apply? If the concept of loop bandwidth is applicable, how do you set it and what are the appropriate values (in relation to symbol rate, maximum frequency offset etc)?

UPDATE

Harris' ML structure

There is still a loop with an adjustable loop bw: The summing block shown is an indication of the additional integration to create a second order loop, although a proportional path should be included as well in order to maintain stability (PI Loop Filter, here with what is suggesting an accumulator you would only have "I"). I suspect the block diagram is just simplified and the summation shown is intended to be the more detailed proprtional-intergral loop filter structure as shown in your subsequent diagram.

Even without the summing block this will still be a loop with adjustable loop BW, just not the best choice: in this case it would be a first order loop that will lock to constant error given a constant frequency offset (and be unconditionally stable if not for parasitic delay elements the implementation will invariably have). If the loop gain is large (set by scaling the value at the input to the DDS), then the error can be minimized- but still exists. Adding the additional proportional-integral loop filter before the DDS input will maintain stability and lock to zero error in this frequency offset condition.

(The DDS, similar to a VCO is directly proportional to frequency- so in units of phase it  forms one integrator in the loop- with the PI loop filter added, we get the second integrator, and with both at DC to form a second order Type 2 Loop).

As I mentioned in the comments, if you do implement this as given, compare results to a fixed scaling (instead of the SNR) and the sign of the signal (instead of tanh) and let us know your results of that versus the optimal implementation shown.


What do the filter coefficients in a digital filter represent? I designed a digital filter using fdatool of matlab and obtained the filter coefficients from the tool.

The problem is that i designed a 4th order filter. This gave me 5 filter values

h[] = {0.1930,0.2035,0.2071,0.2035,0.1930}
x[k] = Discrete time input signal

Now on using the formula

Output = h[k]*x[n-k];

Output represents the final filtered value.Although the results are coming fine, but I am not able to find out how those coefficients are obtained by matlab and how mere multiplication(convolution) gives the final filtered response for any sample.

Any link or explanation will do. I wish to know the complete back-end working of filter coefficient calculation.

Please comment if i am unclear in my doubt somewhere.

Thanks :)

We can try a very short introduction:

Every filter represents a Linear Time Invariant System (LTI)
Every Linear Time Invariant System can be completely described by it's transfer function or it's impulse response. The two can be converted into each other by the Fourier Transform
Filter coefficients are derived from impulse response or transfer function
The exact nature of the filter coefficients depends on the algorithm (there are quite a few of those)
In the case of the simplest algorithm, the direct convolution FIR (Finite Impulse Response) filter, the filter coefficients are simply the impulse response of the LTI system.
In most other algorithms the relationship is much more complicated and text book study is indeed required.
The whole subject of LTI systems, transfer functions, Fourier Transforms, amplitude responses, phase responses etc. is probably another text book worth of stuff


Is it possible to add some kind of noise to a non stationary (having unit root) signals to make them stationary? I am working on time-series signals for which I need to fit a model (like MVAR) that could describe the process satisfactorily. However, my signal has a unit root and is highly non-stationary because of which it doesn't meet the required stationarity criteria for MVAR. With highly non-stationary, I mean it is of differencing order of 7 which makes the interpretation difficult. So, I was wondering if there is any way to make the signals stationary by adding some noise (preferably white) of some specified frequency range.
Thank you!

One processing approach is to use the Two-Sample Variance (also called the Allan Variance), which takes the differencing over intervals of increasing duration from which it can be directly observed over what interval the original signal (unprocessed) can be assumed to be stationary, and for other durations this processing can convert non-stationary signals to be stationary for further statistical evaluation. This is used widely in the Clock world due to the non-stationary characteristics of oscillator phase noise, but is applicable to other signals as well that have similar non-stationary noise processes (1/f noise, drift, etc).
More information on the Allan Variance (or it’s square root the Allan Deviation, ADEV) is here at Wikipedia:
https://en.wikipedia.org/wiki/Allan_variance
Differencing is a high pass comb filtering response, and works in the case of converting non-stationary clock frequency noise since the non-stationary aspects of the noise are more pronounced as you get closer and closer to the carrier. You can also see this with daily stock market prices where the cost vs time is a random walk process (non-stationary) but the daily difference would be a stationary white Gaussian noise process (assuming no memory in the market which is typically valid for prices from one day to the next- even from one minute to the next which is why traders are concerned with  ns time intervals with regards to trading information / communications)
And I have other posts here on DSP.SE that demonstrate the application and utility of the Allan Variance and Allan Deviation beyond the evaluation of Clocks:
see
How to interpret Allan Deviation plot for gyroscope?
and
What determines the accuracy of the phase result in a DFT bin?


Lattice structure realization of an all-pole system Given an IIR system with the system function,

$$ H(z) = \frac{1}{2 + 1.4z^{-1} + 1.8z^{-2}} $$

I need to draw the lattice structure realization.

I have drawn the lattice structure for $\frac{1}{1 + 0.7z^{-1} + 0.9z^{-2}}$, which is shown in the figure below.

Now since $H(z) = \frac{1}{2(1 + 0.7z^{-1} + 0.9z^{-2})} $, should I multiply the signal $x(n)$ by $2$ or $1/2$ before feeding to the system? My intuition is that I should multiply with $\frac 1 2$, but then since $2$ appears as the multiple in the denominator polynomial, I got confused.

Rearrange into a difference equation form and your answer will become clearer!


General questions on the CMA equalizer So I recently implemented a CMA equalizer in MATLAB that uses method of steepest descent to converge to the minimal cost. (Im on an equalizer implementation binge).

My questions are the following:

1) It seems to me that the CMA algorithm is only good for channels where they are relatively 'flat'. In other words, it is not expected to work for channels with deep fades/nulls? Is this generally true?

2) I am using a BPSK signal, and we can see from the first figure here, that after the multipath effects I have a smearing on the complex plane of my BPSK signal - there arent two nice clusters as one would expect. Instead, we see 4 red clusters. My question is how come after the CMA equalizer I still have 4 clusters? (I colored those blue). I suppose it makes some sense because the CMA simply forces the envelopes to be 1, and doesnt 'care' which cluster you are talking about. I have heard that CMA can suffer from local minima problems however. Would this be an example of it? (ie, getting 4 clusters instead of 2 since this is BPSK). If not, what can be done about it?

3) Almost as if answering question 2, I went ahead and changed the constant modulus that I seek to minimize the error against. Instead of choosing 1 (as is supposed to be the case for BPSK), I chose 0.25 to be the modulus. This is the constellation I got:

The problem is that even if this is a 'solution', how does one know apriori what to chose the modulus to be? The reason I consider it a problem is that if I have 4 clusters instead of 2 it makes post-symbol phase/frequency offset estimation/correction more complicated especially when one expects 2 clusters due to BPSK signalling.

(For completeness I have attached the same plots but when I added frequency offsets)

Thanks in advance for any insights you can give into this equalizer!

In general equalizing a channel with deep fades is a problem for all equalizers. CMA equalizer is no exception. Absence of a training signal makes things worse.
This is obviously an example of the CMA stuck in a local minimum. The inital condition of the CMA equalizer was probably in the vicinity of the local minimum. The equalizer filter to which the CMA converged, convolved with the channel does not yield to an impulse-like overall response.
Remember that the cost functions of CMA equalizers are nonlinear. Modifying the modulus do not have a "linear" effect on the cost function. What happened probably was that your initial condition (which is identical to the one in the previous case) of the "modified modulus" CMA equalizer had now fallen into the vicinity of the global minimum. This is the "beauty" of nonlinearity!


Calculate bandwidth of a signal for Nyquist–Shannon sampling theorem I have to calculate the minimum sampling frequency for the Nyquist–Shannon sampling theorem which is Fc > 2*B, where B is the signal bandwidth.
I have this signals: $\text{sinc}^5(t/2 - 4)$ and $\text{sinc}^3(3 - 2t)$

How can I calculate the signal bandwidth to obtain the minimum sampling frequency for the Nyquist–Shannon sampling theorem?

So apply the following principals to get to a solution:

Shifts in the time domain don't affect the bandwidth, so the $-4$ and the $+3$ can be ignored.
Reversal of the time domain axis doesn't affect the bandwidth, so you can ignore the negative sign on the $-2t$.
The $\mathrm{sinc}()$ function has a Fourier transform that is a rectangle function of a particular width.
Scaling the time axis (lengthening or contracting) of a function, scales the frequency axis in the opposite manner (contracting or lengthening) for the Fourier transform.
Multiplication of functions in the time domain, results in the convolution of those two functions' Fourier transforms in the frequency domain.
The convolution of two Fourier transforms in the frequency domain, will result in a Fourier transform whose bandwidth is the sum of the bandwidth of those initial two Fourier Transforms.

That should give you enough to come up with an answer.


What is the physical interpretation of the dB scale on a bode plot and what is a negative dB? I have no physical interpretation of the Bode plot. What does it mean for a bode plot to have negative dB over its entire duration on the log-scale frequency?

Decibels (dB) are used to represent a power ratio with a logarithmic scale. Specifically, a power ratio can be expressed in dB as follows:

$$
R|_{dB} = 10 \log_{10}R = 10 \ \log_{10}\frac{P_1}{P_2}
$$

What does a negative number of dB mean? Manipulate the above equation a bit:

$$
R|_{dB} = 10 \log_{10}R
$$

$$
\frac{R|_{dB}}{10} = \log_{10}R
$$

$$
10^{\frac{R|_{dB}}{10}} = R
$$

If the ratio measured in dB is less than zero, then the exponent on the left hand side will be negative. 10 raised to a negative power results in a number that is in the range $(0, 1)$.

Therefore, if a power ratio measured in dB is less than zero, this implies that, when measured on a linear scale, the ratio is less than one. For a Bode plot, that would mean that the frequency response in question has an amplitude response that is less than unity (for all frequencies, if the Bode plot measured in dB is less than zero everywhere).


How to evaluate expressions without explicitly computing a DFT? I am trying to learn this part but I dont seem to understand.
if I have a finite length sequence, and a n-point DFT of it (in an interval)... is it possible to evaluate expressions of the  n-point DFT without computing them explicitly?

the exercise is:

Let $x[n]$ be finite-length sequence such that: $x[n] = \{2,1,1,0,3,2,0,3,4,6\}$ for $0<=n<=9$
    and $X[k]$ be the 10-point DFT of $x[n]$ for $0\le k\le9$. evaluate the following without explicitly computing the DFT:

$$I) \space  X[0]$$
$$II) \space \sum_{k=0}^9 X[k]$$
$$III) \space \sum_{k=0}^9  |X[k]|^2 $$
$$IV) \space \sum_{k=0}^9 X[k]e^{-j \frac {(4\pi)}{5}k}$$

and then it continues:

let $X[k]$ denote N-point DFT of an N-point sequence $x[n]$ ($n$ is even). Two N/2-point sequence are defined as:

$$g[n] = a_1  x[2n] + a_2  x[2n+1]  \space where \space  0\le n\le(N/2)-1$$
$$h[n] = a_3  x[2n] + a_4  x[2n+1]  \space where \space  0\le n\le(N/2)-1$$

constant $a_i$ satisfies: $a_1  a_4 \neq a_2  a_3$

let $G[k]$ and $H[k]$ denote the N/2-point DFT's of $g[n]$ and $h[n]$ respectively.
Determine N-point DFT $X[k]$ in terms of $G[k]$ and $H[k]$ .

any clues on the solution?

any input is highly appreciated.
thank you

Parts of the answer.

I) is the average value (DC coefficent) of the signal, 2.2

III) is, by Parzeval's theorem, equal to the square sum of the signal, 80


"plot" command for continous time signal in matlab?uses interpolation? We know digital computer cannot deal with continuous time signals directly. But then how MATLAB "plot" command works?

Does that uses interpolation technique to show an apparently continous time graph by connecting the various discrete data points?

Till now, i am able to understand is that MATLAB only handles discrete time data, MATLAB models continous time signals using discrete time signals but with very very small sampling interval,so it appears to be continuous while in reality it is discrete and then connects those discrete points using interpolation

I have also attached a snapshot of wikipedia where a graph of polynomial interpolation is shown

I have also attached a snapshot of chapter4 first page of book signal processing first, where interpolation discussion is mentioned

By default MATLAB uses linear interpolation when creating line plots, which means it simply draws a line from each point to each point, unless there are more points than pixels in which case each point or group of points within each pixel would represent each pixel shown.

Compare the following to see this:

This plots a line from x,y=1,2 to x,y= 5,7

plot([1 5], [2 7])

The above is the same as the following where a line is explicitly declared:

plot([1 5], [2 7], '-')

While in comparison the following will plot a “o” at each sample:

plot([1 5], [2 7], 'o')

Update: To interpolate more smoothly between samples such as polynomial interpolation, you would need to create the interpolated samples and then plot the higher sampled result. The simplest way to do this in MATLAB is to use the resample command where resample(f,M,N) will interpolate the waveform $f$ by $M/N$ samples. Another approach demonstrated below is interpolating using the FFT. Although DSP.SE is not a MATLAB coding site, the following does demonstrate well how to properly do a time domain interpolation using a zero-padded IFFT which is of general signal processing interest. Zero padding in one domain results in interpolation in the other domain, however the technique is to "zero-pad" from the center of the frequency domain waveform with index k from 0 to N-1, such that the high frequency component of the FFT is extended.

N=6;
M=1000;
x= [0:N-1];
y = sin(2*pi*x/N);
fy = fft(y);
padfy = [fy(1:floor(N/2)) zeros(1, M-N) fy(floor(N/2+1:end))];
figure;
plot(N*[1:M]/M, M/N*real(ifft(padfy)));


Simultaneously minimize and maximize two cost functions I want to solve the following optimization problem, in which I have two cost functions $J_1 = \mathbf{w}^T R_1 \mathbf{w}$ and $J_2 = \mathbf{w}^T R_2 \mathbf{w}$, where $R_1$ and $R_2$ are covariance matrices (symmetric PSD).

The problem is to find $\mathbf{w}$ such that $J_1$ is minimized and $J_2$ s maximized constrained on $\mathbf{w}^T \mathbf{w} = 1$.

NB:- This problem arises in signal processing algorithm development and hence posted in this forum.

Thanks for any help.

Purna.

Ideas:

1) Find the minimum eigen vector of $R_1$ and assign this to $w$. This will minimize $ J_1$ , doesn't maximize $J_2$. But when $R_1$ and $R_2$ are positive definite or full rank matrices and computation is an issue, this is a decent solution.

2) Form a new objective $J_1 - J_2$ and minimize this with the given constraint

3) Try and formulate as max SINR problem with $J_1$ as the interference term and $J_2$ as the signal term and the given constraint, you should be able to apply the relaxation $w^Tw < 1$ to convert this to a convex problem. This is readily solved in literature ( simply google (Max SINR) optimization.


inverse discrete fourier transform with plain python I am trying to calculate inverse discrete fourier transform for an array of signals.

I am using the following formula:

$$ x[n] = \tfrac1N \sum\limits_{k=0}^{N-1} X[k] \, e^{j 2 \pi k n/N} $$

And my python code looks as follow.

def IFT(array):
    array = np.asarray(array, dtype=float)
    # array length
    N = array.shape[0]
    # new array of lenght N [0, N-1]
    n = np.arange(N)
    k = n.reshape((N, 1))
    # Calculate the exponential of all elements in the input array.
    M = np.exp(2j * np.pi * k * n / N)
    return 1 / N * np.dot(M, array)

y = 2*np.sin(2*np.pi*f1*t+0.2) + 3*np.cos(2*np.pi*f2*t+0.3) + noise*5
IFT(sc.fft(y))

array([ 3.42732136e+00+0.00000000e+00j,  4.81582993e+00-4.10374772e-16j,
    5.86456023e+00+3.88824673e-16j,  4.07942919e+00-1.48609552e-15j,
    5.25077594e+00+6.70968395e-16j,  6.45823471e+00-1.22977116e-15j,
...
sc.ifft(sc.fft(y))

array([ 3.42732136-8.59756710e-16j,  5.34401724+4.97379915e-17j,
    5.46939384-9.87654403e-16j,  3.31504852-4.17443857e-16j,
    6.52888175+6.39488462e-16j,  6.90804001-3.10862447e-16j,
...

but if I compare array returned from Scipy.ifft with the one returned from my own function IFT, I get completely different array. Am I missing something?

Removing the 2nd term -> array = np.asarray(array, dtype=float) should work


Fourier Series Coefficients Question:

The fourier series coefficients is given as:

$$c_k= \begin{cases}
 1 \qquad & k \ \text{ even} \\
 2 \qquad & k \ \text{ odd} \\
\end{cases}$$

the period of the signal is $T=4$, what is signal $x(t)$?

Attempt:
when i am trying to find the signal by applying the general formula at the end i am getting exponential term now i am finding that hard to convert it into impulse again because i know answer of this signal will be a impulse but i am not getting any idea to convert it into impulse again? Please help with this?

You should use the synthesis equation of an impulse train with period $T$ (which is easy to derive):

  $$x(t)=\sum_{k=-\infty}^{\infty}\delta(t-kT)=\sum_{k=-\infty}^{\infty}\frac{1}{T}e^{jk\frac{2\pi}{T} t}\tag{1}$$
  That is: the Fourier coefficients for all terms is a constant ($\frac{1}{T}$).

Now assume that there are two impulses with different amplitudes $a$ and $b$ per period, or $$x(t)=a\delta(t)+b\delta(t+2),\ -4< t \le0$$ In such case we have

$$\begin{align}
x(t)&=\sum_{k=-\infty}^{\infty}a\cdot\delta(t-kT)+b\cdot\delta(t+\frac{T}{2}-kT)\\&=\sum_{k=-\infty}^{\infty}a\cdot\delta(t-4k)+b\cdot\delta(t+2-4k)\\
&=\sum_{k=-\infty}^{\infty}\frac{a}{4}e^{jk\frac{2\pi}{4} t}+\sum_{k=-\infty}^{\infty}\frac{b}{4}e^{jk\frac{2\pi}{4} (t+2)}\\
&=\sum_{k=-\infty}^{\infty}\frac{1}{4}e^{jk\frac{\pi}{2} t}\left(a+be^{jk\pi}\right)\\
&=\sum_{k=-\infty}^{\infty}\frac{1}{4}e^{jk\frac{\pi}{2} t}\left(a+b(-1)^k\right)\\
&=\begin{cases}
\displaystyle\sum_{k=-\infty}^{\infty}\frac{a+b}{4}e^{jk\frac{\pi}{2} t},& k\text{ even}\\[10pt]
\displaystyle\sum_{k=-\infty}^{\infty}\frac{a-b}{4}e^{jk\frac{\pi}{2} t},& k\text{ odd}
\end{cases}
\end{align}$$
Now refering again to $(1)$ and comparing it with the question, we should have
$$c_k=\begin{cases}
\frac{a+b}{4}=1\\[10pt]
\frac{a-b}{4}=2
\end{cases}$$


find impulse response from step response if we have LTI system and we know unit step response of this system(we haven't original signal)
how we can calculate impulse response?

You can find the impulse response

Let's take the case of  a discrete system.

If $s[n]$ is the unit step response of the system, we can write

$$s[n]= u[n]\ast h[n]$$

where $h[n]$ is the impulse response of the system and $u[n]$ is the unit step function.

Now using commutative property you can write $$s[n]=h[n]\ast u[n]$$

Expanding convolution we get $$s[n] = \sum_{k=-\infty}^{\infty}h[k]u[n-k]$$

which can also be written as

$$s[n]=\sum_{k =-\infty}^{n} h[k]$$

So we can see that unit step response is like an accumulator of all value of impulse response from $-\infty$ to $n$. So now impulse response can be written as the first difference of step response.

$$h[n]=s[n] - s[n-1]$$


How to sketch this multirate system spectrum? In the following system:

The signal x(n) is upsampled by 2 and then filtered. The signal and the filter have the following spectrum:

After upsampling the spectrum of X would look the same except the x axis would be divided by 2. So the limit of the spectrum would be at pi/4. Since the filter has amplitude 2 would the amplitude of X be increased from 1 to 2 or will the amplitude of X remain unchanged?

Finally, for the multiplication of the cosine it would be equivalent to multiply with 2 deltas located at -pi/2 and pi/2. Would the multiplication with this deltas cause a kind of mirroring effect on the spectrum?

Thank you

edit:

I do not see how to obtain Y1 and Y2
complete system

and provided answer:

You are on the right track.

  After upsampling the spectrum of X would look the same except the x axis would be divided by 2. So the limit of the spectrum would be at pi/4.

Yes, that's correct.

  Since the filter has amplitude 2 would the amplitude of X be increased from 1 to 2 or will the amplitude of X remain unchanged?

Recall that in the frequency domain, the application of a linear system transfer function $H_B(\Omega)$ is done via straight multiplication of the input signal's spectrum. So yes, since the filter has a magnitude of 2, the amplitude of its output is doubled.

  Finally, for the multiplication of the cosine it would be equivalent to multiply with 2 deltas located at -pi/2 and pi/2. Would the multiplication with this deltas cause a kind of mirroring effect on the spectrum?

Here's where I will make a small correction: multiplication by a cosine in the time domain yields convolution with a pair of Dirac impulse functionss in the frequency domain. Convolution with a Dirac impulse in the frequency domain is equivalent to a shift along the frequency axis. So, you end up with 2 copies of the filter's output spectrum, centered at $\omega = \pm \frac{\pi}{2}$.


Use of Tukey window and its effect on a graph I am trying to minimise the noise in a spectrum of a variable. I am trying to to divide the time history in equal parts, to do the FFT of this time history parts and to do an average between these fft.
Doing this procedure I have find that the trend, at the high frequencies seems to be a little bit different from the one without doing this procedure. Following an other suggestion of my advisor I have multiplied every part of the time history by a Tukey window, using the command with matlab
A= tukeywin(L)
Using the Tukey window, I have find that the trends seems more similar with the first one. Can someone explain me why? I have tried to search the theory of this window, but I always find the definition of the filter, but I cannot understand which could be the effect of this window.
A other fact that I have noticed is the fact that obviously using the tukey window the energy content is reduced, are there technique to reduce the noise avoiding this reduction?
I just need a suggestion to start, it is difficult to understand which is the correct technique to reduce the noise.

The technique that you have described looks similar to the Welch's method for estimation of the power spectral density. Though, it seems that there is an error in your implementation: you should take the squared absolute value of the FFT, instead of the complex value. Such an error might affect your results significantly as the random phase might bring you to a zero mean value of the spectrum. Anyway, most of the languages used for signal processing such as MATLAB and Python contains a ready to use implementation o the Welch's method, so you do not have to mess with the details.
Suggestion: since the variation of the spectrum is due to the user of several windows (variation reduction proportional to N), you might want to overlap the windows. The typical overlap for the common Hann window is 50%.

Window type should not affect the spectrum in a way that you have described. Multiplication in time is equivalent to convolution in the frequency. Therefore, the windowing effect should mostly affect your spectrum leakage. Is your signal stationary? this is a fundamental assumption for using such a technique.

Actually, as long as I know, the signal that you have described does not contain information in the time domain but only in the frequency domain. In other words, there should be no periodic components with defined phase. Therefore, it is usually processed with a rectangular window (and then there is no need for overlap).

Using any window reduces the energy as the window's value is smaller than 1. This could be compensated by dividing the result with the window norm/mean value (maybe with the power of two or some other coefficient - depends on the type of the spectrum you want to estimate). Anyway, you don't have to worry about it if you use the implementation by MATLAB which wrappers it all for you.


Comparing band-pass filters and convolutional kernels? I am currently working on extracting features for speech recognition purposes.
I wanted to try an approach similar to MFCC features in which the center frequency of each band pass filters are placed a distance of x mel displaced of each other.

But instead of looking at the audio in 1d I see it in 2d using the spectrogram as representation, and I was thinking of using a CNN instead of using band pass filters , with different weight for sets of frequencies, as they are interpreted differently, and some more distinct than others.

Elaborate on some are more distinct than others:

Consider the normal mel scale.  Frequency between 0 - 1000 are interpreted linearly increasing, and above as logarithmic increasing. The way I interpret this is that the ear is more susceptible to differences in the lower range (<1000), and not so much in the higher range (>1000). So in my head it makes sense to place the weight as higher in the lower range, and decreasing the higher range.

But what about the kernel size and the center frequency?.. As I enter the logarithmic range, will the band pass filters used in the MFCC become wider, which in this case should result in larger kernels, but what about the center frequency?.. Where or how should i place a kernel such that a term such as center frequency would make sense related to kernels?... is it even possible, or am I just being silly?

CNNs sort of irrelevant here because they are applicable when you have some invariant transform - shift or rotation. The idea of CNN is that weights are shared for different data points, not that weights are different. You are thinking about simple DNN.

CNNs are also applied in speech recognition, but for a bit different task. For example, you can consider multiple frames jointly and tie filter shape for different frequencies. This way system learns Gabor-like filters. See Robust CNN-based Speech Recognition With Gabor Filter Kernels.

You can try to train your own kernels from raw data, they will eventually resemble bandpass filters. See Acoustic Modeling with Deep Neural Networks Using Raw Time Signal for LVCSR. This gives a marginal improvement.

You can also start from default bandpass filters as initial values and see if network learns an improved version.

You can also restrict network architecture to consider only the region around central frequency and just learn the weights in that region. Such network will unlikely improve anything.


how is this method of proving shift invariance correct I got this answer from here. But i don't know how this method is correct. Until now to prove shift invariance what i did was :
step 1: delay output and write it in terms of input.
step 2: delay input and write the modified output
step 3: if output from step 1 and 2 are same then system is shift invariant.

what this guy has done is "different". the oppenheim book and proakis manokalis book don't follow this method. which book follows this method?

the method.

  The test for this is
  $$x_{1}[n]=x[n-n_{0}] $$
  $$y_{1}[n]=y[n-n_{0}] $$

  So let's first express in terms of just $y[n]$.
  $$ y[n]=\frac{x[n-3]-y[n-3]-y[n-1]}{x[n]} $$
  Next, we go through the test.
  $$y[n-n_{0}]=\frac{x[n-3-n_{0}]-y[n-3-n_{0}]-y[n-1-n_{0}]}{x[n-n_{0}]} $$
  Substituting $y_{1}[n] $
  $$y_{1}[n]=\frac{x[n-3-n_{0}]-y[n-3-n_{0}]-y[n-1-n_{0}]}{x[n-n_{0}]} $$
  $$y_{1}[n]=\frac{x[n-3-n_{0}]-y_{1}[n-3]-y_{1}[n-1]}{x[n-n_{0}]} $$
  And finally substituting $x_{1}[n] $
  $$y_{1}[n]=\frac{x_{1}[n-3]-y_{1}[n-3]-y_{1}[n-1]}{x_{1}[n]} $$
  Because the shifted sequence has the exact same relationship, it is said to be time-invariant. Most systems that don't alter $n$ or $t$ meet this. Be wary of anything multiplying or otherwise messing with the vector arguments beyond simple delays.

The so called example system is nonlinear and recursive and those books generally exclude the treatment of nonlinear & recursive systems. Hence you may not find example solutions.

The method relies on the fact that the relationship between the output $y[n]$ and input $x[n]$ remains the same between $y[n-d]$ and $x[n-d]$.

As an absurd example, think of the following. Consider the equation of a circle at origin: $x^2 + y^2 = r^2$. Now, shift this circle to the point $(a,b)$ and re-write its equation: $(x-a)^2 + (y-b)^2 = r^2$. Finally, redefine the variables $x' = x-a$ and $y' = y-b$ and rewrite the equation: $(x')^2 + (y')^2 = r^2$.

It's seen that the circle is shift invariant (not in the sense of the signal processing context though)


Underdetermined deconvolution of windowed output Consider a discrete 'blurred' output $h[t]$ given by the convolution of filter $f[t]$ and signal $g[t]$. This question considers recovering $g[t]$ from a window (subset) of $h[t]$. This causes the problem to be under-determined.

$$h[t] = (f * g)[t]$$

Original signal $g[t]$:

$$g[t] = [1,5,3,5,7,7]$$

Moving average filter $f$ of window size 2:

$$f[t] = [0.5,0.5,0,0,0,0,0]$$

The full 'blurred' output $h$ is given by:

$$h[t] = [0.5,3,4,4,6,7,3.5]$$

Given $h$ the problem of recovering $g$ is well-posed and can be recovered simply using Fourier transforms:

$$g[t] = \mathcal{F}^{-1}\left\{\frac{\mathcal{F}\left\{h[t]\right\}}{\mathcal{F}\left\{f[t]\right\}}\right\}$$

However in practice the recovery needs to be made from a windowed output $h'[t]$:

$$h[t]' = [3,4,4,6,7]$$

The problem of recovering $g[t]$ from $h'[t]$ is ill-conditioned since there are an infinite number of solutions. One approach I have seen suggested is to take $g[1] = h'[1]$. In this case the recovered profile $g'[t]$ is given by:

$$g'[t] = [3,3,5,3,9,5]$$

My understanding is that most approaches are based on regularisation however most sources seem to focus on ill-conditioning resulting from additive noise rather than ill-conditioning resulting from a windowed output. What are the approaches typically used for this problems? Can you apply it to this toy example?

Using more conventional notation, let $x_k$ and $y_k$ denote the $k$-th input and $k$-th output, respectively.

$$\begin{array}{rl} y_1 &= 3 = \frac 12 x_1 + \frac 12 x_0\\ y_2 &= 4 = \frac 12 x_2 + \frac 12 x_1\\ y_3 &= 4 = \frac 12 x_3 + \frac 12 x_2\\ y_4 &= 6 = \frac 12 x_4 + \frac 12 x_3\\ y_5 &= 7 = \frac 12 x_5 + \frac 12 x_4\end{array}$$

We have an underdetermined system of $5$ linear equations in $6$ unknowns. Let $x_0$ be a parameter.

$$\begin{bmatrix} \frac{1}{2} & 0 & 0 & 0 & 0\\ \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0\\ 0 & \frac{1}{2} & \frac{1}{2} & 0 & 0\\ 0 & 0 & \frac{1}{2} & \frac{1}{2} & 0\\ 0 & 0 & 0 & \frac{1}{2} & \frac{1}{2}\end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3\\ x_4\\ x_5\end{bmatrix} = \begin{bmatrix} 3\\ 4\\ 4\\ 6\\ 7\end{bmatrix} - \frac 12 x_0 \begin{bmatrix} 1\\ 0\\ 0\\ 0\\ 0\end{bmatrix}$$

Using SymPy, we can solve the linear system

>>> from sympy import *
>>> from fractions import Fraction
>>> A = Fraction(1,2) * Matrix( [[1,0,0,0,0],
                                 [1,1,0,0,0],
                                 [0,1,1,0,0],
                                 [0,0,1,1,0],
                                 [0,0,0,1,1]] )
>>> y = Matrix([3,4,4,6,7])
>>> x0 = Symbol('x0')
>>> A**-1 * (y - Fraction(1,2)*x0*Matrix([1,0,0,0,0]))
Matrix([
[-x0 + 6],
[ x0 + 2],
[-x0 + 6],
[ x0 + 6],
[-x0 + 8]])

Thus, the solution set is a line parametrized as follows

$$\mathrm x \in \left\{ \begin{bmatrix} 6\\ 2\\ 6\\ 6\\ 8\end{bmatrix} + x_0 \begin{bmatrix} -1\\ \,\,\,\, 1\\ -1\\ \,\,\,\, 1\\ -1\end{bmatrix} : x_0 \in \mathbb R \right\}$$

Note that if we choose $x_0 = 1$ then we recover the original input vector. You chose $x_0 = 3$ instead.


What should be considered when selecting a windowing function when smoothing a time series? If one wants to smooth a time series using a window function such as Hanning, Hamming, Blackman etc., what are the considerations for favouring any one window over another?

The two primary factors that describe a window function are:

Width of the main lobe (i.e., at what frequency bin is the power half that of the maximum response)
Attenuation of the side lobes (i.e., how far away down are the side lobes from the mainlobe). This tells you about the spectral leakage in the window.

Another not so frequently considered factor is the rate of attenuation of the sidelobes, i.e., how fast do the sidelobes die down.

Here's a quick comparison for four well known window functions: Rectangular, Blackman, Blackman-Harris and Hamming. The curves below are 2048-point FFTs of 64-point windows.

You can see that the rectangular function has a very narrow main lobe, but the side lobes are quite high, at ~13 dB. Other filters have significantly fatter main lobes, but fare much better in the side lobe suppression. In the end, it's all a trade-off. You can't have both, you have to pick one.

So that said, your choice of window function is highly dependent on your specific needs. For instance, if you're trying to separate/identify two signals that are fairly close in frequency, but similar in strength, then you should choose the rectangular, because it will give you the best resolution.

On the other hand, if you're trying to do the same with two different strength signals with differing frequencies, you can easily see how energy from one can leak in through the high sidelobes. In this case, you wouldn't mind one of the fatter main lobes and would trade a slight loss in resolution to be able to estimate their powers more accurately.

In seismic and geophysics, it is common to use Slepian windows (or discrete prolate spheroidal wavefunctions, which are the eigenfunctions of a sinc kernel) to maximize the energy concentrated in the main lobe.


Are resolution increase and noise reduction from oversampling mutually exclusive? Oversampling a signal means sampling it with a significantly higher sampling frequency than the Nyquist rate. As far as I know, there are three advantages:

Easier design of anti alias filter
Increase in resolution in the presence of a white noise signal (I was told to assume the quantization error can under the right conditions be taken as this noise). While the uncorrelated noise adds up destructively (amplitude raises by $\sqrt N$), the signals amplitude increases by a factor of $N$. This means the SNR increases by $\sqrt N$, acting as an increase in resolution.
Decrease in noise. The quantization noise power is independent of the sample rate and only depends on the amount of bits. So the same noise power is distributed over a larger frequency band, resulting in less noise power per frequency. The a low pass filter can be used to remove the high frequency part of the signal.

If I average over multiple values I am back to the lower sample rate I started with, so I don't think the noise would decrease as described in 3. This leads me to think 2. and 3. are exclusive.
My question is, are 2. and 3. exclusive or can you benefit from both when oversampling?

If you just average multiple samples, you are using a very low quality anti-aliasing low-pass filter which can add a ton of noise. If you then re-quantize after averaging, you add still more noise.

If you don't (re)quantize after a high quality anti-aliasing low-pass filter (but keep a ton more bit of resolution), then noise could still decrease as per (3).  And if you FFT before your resampling to a lower sample rate, you still get the higher bin resolution of a longer FFT.  So (2) and (3) need not be mutually exclusive.


Get data ready for IFFT I am having troubles applying the IFFT function to some data I have. I already had a look here Complex conjugate and IFFT and here What assumptions should be used to invert spectrum into time domain data? but it still does not work.

This is the code:

% amplitude: data related to the signal amplitude
% phase: data related to the signal phase
% freq: frequencies (from 5 to 100 Hz)

response = complex(amplitude,phase);

NFFT = length(response); % odd number
symmetric = zeros(NFFT,1);

for i = 1:NFFT
symmetric(i) = conj(response(mod(NFFT-i+1,NFFT)+1));
end

response_tot = [response; symmetric];

% Plot simmetry
response_tot = fftshift(response_tot);
f_plot = 100*(-NFFT:(NFFT-1))/NFFT;
figure()
plot(f_plot,response_tot)

% Time signal
time_sig = ifft(ifftshift(response_tot));

fs = 2*max_f; % twice the max frequency in the signal
t_plot = 0:1/fs:NFFT/fs; % time vector to plot the time response

figure()
plot(t_plot,time_sig(1:NFFT+1))

Is this the right way to prepare the data for the ifft?

Because the output does not look right to me.

The first component of the FFT is the DC component. The other NFFT-1 should be symmetric.
If NFFT is even, that means that

response(NFFT/2-k)=conj(response(NFFT/2+k))

for $k=1,\ldots,NFFT/2-1$ and response(NFFT/2) is arbitrary.
If NFFT is odd then

response((NFFT+1)/2-k)=conj(response((NFFT+1)/2+k+1))

for $k=0,\ldots,(NFFT-1)/2-1$.

However, MATLAB allows you to ignore all that by using the 'symmetric' argument of ifft:

response_tot = [response; zeros(NFFT-1,1)]; % assuming NFFT is odd
time_sig = ifft(response_tot,'symmetric');


Phase Spectrum in Fourier Transform Probably a simple question to answer. I have been trying to use the fft function in MATLAB and have been succesful to get the amplitude spectrum of a known signal. But when working out the phase spectrum I am unsure if its the correct output.
Can anyone explain what the phase spectrum in fft is used for? I have researched other places but have not a satisfactory answer.
If for example I had a pure sine wave what would the phase spectrum show and how will I know if it is correct.
I can give example of my MATLAB code if necessary but if someone can give a basic explanation using MATLAB it will be much appreciated.

I think that the plot you display does not show the correct answer, what makes me think that is that you inject in your fft a pure sinewave with an initial phase equal to 0. Therefore your plot should display a 0 phase shift for every frequency bins.
For exemple, below you can see the phase plot of a pure sine wave with an initial phase shift of 30 degrees.

The phase spectrum will be here to tell you in which way your signals are delayed for each frequency components.


The order of processing the prerecorded samples (DC removal, Hilbert transform, frequency shifting) I have a samples prerecorded by a radio device. Samples are not in baseband, contain large DC component and include only one channel (written in int32 format). My goal is to move it to the baseband, remove the DC spike and perform Hilbert transform to make it complex for further processing. My question is, what is the correct order of performing mentioned computations? Or it doesn't matter? I assume that first I need to remove DC by calculating mean (DC component) and correct each sample. Then go with Hilbert transform to get complex representation of the signal. After that I need to multiply the signal by a complex waveform to move it to the baseband. Can anyone confirm that my way of thinking is right or correct me if I am wrong.

Order of operations can be changed but will change where the processing is applied and that may be simplified in certain conditions which will be very obvious once the actual signal is processed if the goals of each step are understood. For example, if the OP means by DC a non-zero carrier (which will appear at DC after moving the signal completely to baseband), this is much easier to remove after moving to baseband by simply subtracting the mean. However in actual applications, the carrier may not initially be perfectly estimated, in which case a tone close to DC will occur and subtracting the mean will be ineffective (but the tone can be used to continue to refine the carrier estimate).  If the OP meant the waveform itself as received prior to any processing has a large DC offset, this too can be removed prior to any processing with a simple subtraction of the mean (since the samples are pre-recorded this would be simplest) assuming the DC offset is not desired for some reason, otherwise the tone it becomes will be filtered out regardless in the subsequent processing below. (There can of course be advantages to removing the DC signal first in terms of dynamic range available but that would all be clear in the specicic procssing done and precision of the operations.
No Hilbert transform is necessarily needed, if the frequency initially is at a carrier sufficiently higher than its bandwidth, simply frequency translate the real signal to baseband using using $e^{-j\omega_c t}$ which means feed the signal into to multipliers and multiply one with $\cos(\omega_c t)$ and multiply the other with $\sin(\omega_c t)$ to get a complex I and Q output and filter the remaining higher frequency component that will reside at $e^{-j 2\omega_c t}$ (this is a typical Digital Downconverter (DDC) architecture).
For further details of this approach, please see:
https://dsp.stackexchange.com/a/51280/21048


unsigned audio formats, does DC offset matter? I'm a little confused about DC offset in audio files. If I have a signal that's not using the entire dynamic range, does the DC offset matter at all, or does it need to be centered around some zero value? For say unsigned 8-bit audio does that mean i should add my signal to 127, or can I just put it anywhere as long as it doesn't clip and that will be filtered out somewhere down the line?

Consider a speaker which wobbles in and out as its driven by a fluctuating audio curve signal ... when the given audio is silent the speaker driver is at rest, not moving ...  that is if the input audio is properly normalized such that it has no DC offset ...  speakers are engineered to deliver optimum fidelity when given zero DC offset audio curve

Now when the speaker is driven by audio which does contain a DC offset the speaker driver is no longer fluctuating in and out centered at its silent resting driver code position ... so when audio is silent yet with a DC offset I would say the speaker driver is held at a physical offset away from its resting position ... at the very least this will reduce the efficiency of the speakers so perceived volume will be lower and I dare say with reduced fidelity

UPDATE googling :  if audio has a DC component offset will silence make speaker cone stay off its resting position ... wisdom is its bad for speakers  https://electronics.stackexchange.com/questions/233230/why-are-dc-signals-bad-for-loud-speakers


Why are recursives methods useful for FIR filter design? As I understand it, FIR filtering is a linear process. That mean for me that the whole filtering process will have a fully predictable behavior.
So, could someone explain why a universal deterministic and optimal filter design method  to obtain the desired response doesn't exist and why recursive approach seems to be of great interest in practice ?

EDIT: just to clarify, my question is about an "optimal" method, that is one which compute the shortest filter response while staying the closest to the design criteria.

Sure, FIR filters are "simple" linear objects. But a typical FIR filter design task looks like this:

  Find a set of numbers (the filter coefficients) such that the
  magnitude of their Discrete Fourier Transform is as close as possible to a
  target response. "As close as possible" being defined as "by
  minimizing the $L_\infty$ norm of the difference between ideal and actual
  response over some frequency ranges".

This is a non-linear, sometimes non-convex optimization problem and there is no universal, trivial procedure to solve those.

We have to use procedures like the Parks-McClellan algorithm because the real-world engineering constraints we want to impose when designing FIR filters do not translate into simple linear or quadratic constraints. For example, a least-square error criterion (L2 norm often leads to more tractable procedures than other norms...) would not be an acceptable choice for engineering applications, because it could allow solutions that are close to the required response on average, but with an outlier value.


Removal of blocking artifacts in images due to compression Having compressed an image,how can the compression noise or blocking artifacts be suppressed to improve the quality of the image before further processing. Thanks

If you talk about JPEG-like blocking artifacts, an interesting appproach consists in the re-applications of the same JPEG algorithm on shifted versions of the image, see Enhancement of JPEG-Compressed Images by
Re-application of JPEG, by Aria Nosratinia, 1999. It has a software implementation. It implement a sort of shift-invariant processing. He did the same with the JPEG-2000 coder, in Post-Processing of JPEG-2000 Images to Remove Compression Artifacts. You can find an overview of these artifacts here.

Improved methods include sparsifying priors like total variation, see for instance Adapted Total Variation for Artifact Free Decompression of JPEG Images. Recent references include Reducing Artifacts in JPEG Decompression Via a Learned Dictionary and Understanding, Optimising, and Extending Data Compression with Anisotropic Diffusion.


Homodyne and heterodyne detection Very basic question here. I searched it on the net but cudn't get very convincing answers.
Can anyone explain me Homodyne and Heterodyne detections in rf domain. Mathematically is much appreciated.

I'll assume AM DSB modulation for simplicity, but the same idea applies to quadrature modulation. I'll assume that the received signal has been properly bandpass filtered. I'll also ignore noise and carrier synchronization.

So, the received signal is $$r(t) = s(t)\cos(2\pi f_c t).$$ We want to recover the baseband signal $s(t)$ from $r(t)$. The bandwidth of $s(t)$ is $B$ and the bandwidth of $r(t)$ is $W=2B$. The signal $r(t)$ is centered on $f_c$; we would like the receiver to accept $f_c$ as a parameter set by the user (we say that the receiver should be tunable).

A homodyne receiver calculates $$\hat{s}(t) = \text{LPF}\lbrace r(t) \cos(2\pi f_c t) \rbrace.$$ While mathematically this is super simple, in practice this receiver is extremely difficult to implement. This was especially true in the early age of radio, where the entire system was analog and the theory was not yet well understood.

The heterodyne receiver (invented by the unsung, tragic genius Armstrong, as pointed out by Stanley) recovers $s(t)$ in two (or more) stages. First it calculates
$$s_{IF}(t) = \text{BPF}\lbrace r(t) \cos(2\pi f_1 t) \rbrace,$$ where $f_1$ is selected so that $s_{IF}(t)$ is centered at frequency $f_{IF}$.  In a second stage, the receivers calculates $$\hat{s}(t) = \text{LPF}\lbrace s_{IF}(t) \cos(2\pi f_{LF} t) \rbrace.$$ The point of this receiver is that the second stage is independent of the chosen value of $f_c$. This single idea simplifies the implementation so much that it is actually the reason we've had radio at all for the past 70-80 years.

If you're interested in this subject, there are two books I can't recommend enough:

"The science of radio", by Paul Nahin.
"Man of high fidelity", a biography of Armstrong by Lawrence Lessig.

On the technical side, I like the presentation of different modern receiver architectures in

"RF system design of transceivers for wireless communications", by Qizheng Gu.


Fourier Transform Signals - Time Transformations I was going over some review problems and came across an interesting one.

Using the techniques of (linearity, time shifting, and time scaling) what are some approaches I could use to turn the fourier transform from the left graph to the right?

X(ω) = 4* (sin w)^2 / w^2

My understanding so far:

Using the time shift property for fourier transforms I know the we can use
x(t -1) to move the graph to start at t=1
We can use x(t/3) to extend or "stretch" the base to the to t=7

What is eluding me is how to form that "middle" or clipped top of the z(t) graph.
Any insight?

HINT:

Note that $x(t)+x(t-1)$ has the required shape. Now it's only a matter of shifting, time scaling, and amplitude scaling to arrive at the exact expression for $z(t)$.


Coherence Calculation in Sparse Sensing i have an image I of size 32*32. I perform the DCT of this image using the matlab function DCT2(I). I get a sparse representation of my image which is again a 32*32 image. I construct a circulant matrix A of size 200*1024, i.e i need only 200 samples out of the total 1024. Mathematically spaeaking, what i did is
                      take my original signal x, sparsify using dct2 call it B, then compress is using A. So i get measurements y=A.B.x . Now how do i calculate the incoherence measure from A and B. Can somebody help.

Update

I would like to rephrase the original problem with slightly different notation and claim that there is no need to measure the coherence of the product of the circulant matrix A and the DCT matrix B as asked by the OP.

In the compressive sensing setting, the coherence of a sensing matrix $A$ is useful when we are using it for constructing measurements $y = A \alpha$ where $\alpha$ is a sparse vector. Then, the coherence properties of $A$ can be used to give an estimate of number of measurements required in $y$ so that $\alpha$ can be recovered stably from $y$. $\alpha$ is actually a sparse representation of some signal $x$ in some orthonormal basis $\Psi$ which is represented by the equation $x = \Psi \alpha$. In this particular case, $x$ is the image 32x32, $\Psi$ is the IDCT2 transform, $\alpha$ is the sparse representation of $x$ obtained by DCT2 transformation which is the equation $\alpha = \Psi^T x$. Thus, the matrix $B$ from OP is $\Psi^T$. So the relationship between the measurements $y$ and the original signal $x$ is
$y = A \Psi^T x$. The sparse recovery of $x$ from $y$ would happen in two steps.

Construct $\alpha$ from $y$ by using some algorithm $\Delta$. We represent this algorithmic recovery step as:  $\alpha = \Delta(A, y)$. $\Delta$ could be basis pursuit or OMP or anything else.
Construct $x$ from $\alpha$ by using the equation $x = \Psi \alpha$.

In this overall recovery process, the matrix $A B = A \Psi^T$ plays no role. Hence, there is no need to measure the coherence of $AB$.

The old answer was incorrect. The coherence of a matrix $A$ doesn't change if it is pre-multiplied by an orthogonal matrix $B$ as $BA$. Coherence of a matrix can change in arbitrary ways if it is post-multiplied by an orthogonal matrix as $AB$. I am not aware of any particular results of relating the coherence of $A$ with $AB$.

Below is old answer.

You need to compute the coherence of the matrix A. B is an orthonormal matrix. So it doesn't impact coherence of the combined matrix AB .  For computing the coherence of A, normalize its columns and then take pairwise inner products of the columns. The coherence is the highest inner product in magnitude of the pairs of normalized columns of A.

For more information about coherence, see my tutorial notes here .

You may also want to learn about Babel function which is a generalization of coherence.


Cyclic prefix length in baseband OFDM with ZF equalizer My question is about baseband OFDM. I try to plot the BER according to the SNR. I noticed that changing cyclic prefix length does not affect the performances. I.e., following the code below, gi = n_sp/1024 or n_sp/4 or even 0, gives the same performance (always the same BER curve). Is there a mistake somewhere?

NB: At the receiver I use ZF-equalizer. The full Matlab code is:

clear
close all
clc

n_subc = 2^11;                          % # of subcarriers
n_ofdm_sym = 10;                        % # of OFDM symbols
n_data_sym = n_subc*n_ofdm_sym;         % # of data symbols to transmit
gi = n_subc/4;                          % CP length

dqpskmod  =  comm.DQPSKModulator;
dqpskdemod = comm.DQPSKDemodulator;

h = [.7*exp(1i*pi/4) .4*exp(1i*pi/3) .2*exp(1i*pi/6) .1*exp(1i*pi/5)]; % 4-taps channel

SNR = -10:2:30;

for i = 1:length(SNR)
    for j = 1:100                                % Monte Carlo iterations

    d_symb = randi([0 3],1,n_data_sym);
    x_mod = step(dqpskmod, d_symb').';           % DQPSK modulated symbols

    x_mod_s2p = reshape(x_mod,n_subc,[]);        % serial to parallel

    x_ifft = ifft(x_mod_s2p,n_subc);

    x_cp = [x_ifft(end-gi+1:end,:); x_ifft];     % add CP

    x_ifft_p2s = reshape(x_cp,1,[]);             % parallel to serial

    % --------------------> Channel mixing and adding noise

    y = conv(x_ifft_p2s,h,'same');                          % Linear convolution
    % y = ifft(fft(h,length(x_ifft_p2s)).*fft(x_ifft_p2s)); % OR circular convolution
    y = awgn(y,SNR(i),'measured');                          % add noise

    % --------------------> Receiver

    y_s2p = reshape(y,n_subc+gi,[]);             % S2P

    y_cp = y_s2p(gi+1:end,:);                    % Removing CP

    y_fft = fft(y_cp,n_subc);                    % FFT

    y_p2s = reshape(y_fft,1,[]);                 % P2S

    pilotX = x_mod(1:4:end);                     % Pilot extraction from the DQPSK transmitted data
    pilotY = y_p2s(1:4:end);                     % Pilot extraction from the DQPSK received data
    H_est = pilotY./pilotX;                      % Channel estimation
    H_est_interp = interp1(1:4:n_data_sym,H_est,1:n_data_sym); % Linear interpolation
    H_est_interp(end) = H_est_interp(end-1);     % To avoid the last NaN

    X_est = y_p2s./H_est_interp;                 % ZF equalizer

    y_demod = step(dqpskdemod,X_est.')';         % DQPSK demapping

    d_symb(1:4:end)  = [];                       % Pilot suppression from the original data
    y_demod(1:4:end) = [];                       % Pilot suppression from the received data

    [~,ber(j)] = biterr(d_symb,y_demod);

    end
    BER(i)=mean(ber)
end

figure
semilogy(SNR,BER)

The figure below gives the results for gi = n_sp/1024, n_sp/4, and 0:

Additional question: why BER saturates for SNR above 22dB? There is no saturation for QAM modulation, is it specific to the DQPSK?

The channel which you have created is having 4-Taps and all taps are one after the other, meaning roughly there are 4 multi-paths and they are very close to each other. How close depends on what is the Sub-Carrier Spacing you would have assumed. Anyway, the point is only 3 samples of Cyclic Prefix would be enough and even with 0 or 2 samples of CP will not cause any noticeable degradation. The degradation is mostly due to AWGN which you are adding.

If you want to see the impact of 4-Tap channel, then I will suggest you to have larger delays multiplied to the magnitude of the four taps. One way to do it is as follows:

h = [.7*exp(1i*pi/4) zeros(1,50) .4*exp(1i*pi/3) zeros(1,75) .2*exp(1i*pi/6) zeros(1,100) .1*exp(1i*pi/5)]; % 4-taps channel

This will create four paths channel with 2nd, 3rd and 4rth path coming at 50, 126 and 227 samples respectively and hence the 2nd, 3rd and 4rth copy of one OFDM symbol will interfere with next OFDM symbol more.

With your current h assignment, all you will have is intra-symbol interference which can be delt with no CP since the channel is just 4-taps with no delays in between.

Also, consider removing AWGN and varying the amount of zeros I have added in between 4-Taps of h, to actually notice the impact of multipath inter-symbol interference and advantage of 512 length CP over 0/2 length CP.


What is the purpose of dividing a signal into subbands when compressing an audio file? I must implement my own very simple MUSICAM algorithm (basic filtering + Huffman coding) to compress an audio file. But I'm having trouble understanding the interest of dividing a signal into 32 frequency bands before filtering the frequencies. I've read stuff about Nyquist rule but it's not clear for me since I don't have a huge background in signal processing.

Could anyone explain why not just filtering the raw data without dividing into subbands ?

Here is the main reference I use in understanding the concept of MUSICAM.

We don't hear all frequencies equally well. Therefore, if we split the frequencies into subbands, we can give more bits to bands which we hear well, and less bits to bands we don't hear that well (which is determined by things like doing listening tests). Its usually easier easier to control the error of the compressed signal on a subband level than on the whole signal level.

This will typically give us less distortion for a given data rate than trying to encode all frequencies with equal number of bits.

An analogous experiment you can do is with (natural) images: Take the DFT or DCT of an image, and throw out the higher frequencies (i.e. high in the x,y directions). You still get a pretty decent approximation if you keep only a few low frequencies. So, instead of throwing them out, imagine you quantized them -- you'd want to give less bits to the higher frequencies than the lower frequencies, since you're not getting much perceptual gain by throwing lots of bits in the higher frequencies, but you are getting more in the lower frequencies. (This underlies JPEG; See something like Gonzalez and Woods' Digital Image Processing for details on how exactly JPEG works).


How to delay an audio signal by 1 microsecond? I have an audio signal sampled at 44100Hz. I want to delay this signal by 1 microsecond (and maybe even less) for a steganographic purpose. Now if I delay it one sample, the corresponding time delay will be $\frac{1}{44100} = 2.2 \times 10^{-5}\,\text{s}$. But the delay I need is even lesser than that. How to achieve it?

It's called a "fractional delay".

A really comprehensive analysis on how to do this can be found here: https://ieeexplore.ieee.org/document/482137/ and the design tools are still available here http://legacy.spa.aalto.fi/software/fdtools/


Fourier transform and impulse function $\delta(\omega)$ Why does impulse function $\delta(\omega)$ keep occurring in the Fourier transform expression of standard functions like $\sin(t)$, $\cos(t)$, constant function, unit step $u(t)$ etc? (can someone intuitively explain what it means)

The Dirac delta impulse $\delta(\omega-\omega_0)$ represents a spectral line at frequency $\omega_0$, since it is zero everywhere except for $\omega=\omega_0$. So any function with spectral lines, such as a sinusoid, or a DC signal (which has a spectral line at frequency $\omega_0=0$) has a Fourier transform which contains Dirac delta impulses. This is also true for any function with a DC component such as the unit step $u(t)$. Note that also periodic functions have Fourier transforms consisting of Dirac impulses because they can be represented as a sum of sinusoids (Fourier series).


Convert sine wave signal to square wave signal I have a signal containing a variation of digital ASK type on-off keyed sine wave among which some stochastic irrelevant harmonic waves occur, like you can see on this figure:

In its current form each pulse duration is the same, but it's not always the case.

I'd like to calculate each pulse duration. To do that, I was thinking of converting this signal to a square wave signal like this:

I'm not sure of the different steps needed for this conversion but I have some leads:

Remove the noise => Low pass filter I presume?
To have the signal between 0 and the amplitude, I multiply the signal by itself
I apply a magic filter to convert the signal to square signal.

So my questions are:

Do I need a low pass filter to remove the noise?
What kind of filter do I need to convert the signal to square signal? Is it possible?

For information, I'm using Python and Numpy.

If there is any other method to calculate the duration of each pulse in the sequence, I'm open to all ideas.

Thanks!

Edit: Here's my signal:

   freq = 1000000
   fs = int(2*freq) # sample rate
   x = np.linspace(0, 10, fs)
   noise1 = 0.01*np.sin(2*np.pi * 1*x) + 0.05*np.sin(2*np.pi * 1.5*x) +     0.04*np.sin(2*np.pi * 0.2*x)
   noise2 = 0.1*np.sin(2*np.pi * 0.1*x) + 0.2*np.sin(2*np.pi * 1.8*x) + 0.1*np.sin(2*np.pi * 0.4*x)

  realSignal = 2*np.sin(2*np.pi * freq*x)
  sig =  noise1

  sig = np.append(sig, [realSignal, noise2, realSignal, noise1*noise2])

What you can do is process this signal with a windowed RMS function (I'm writing in matlab code, but it can easily be translated in python) :
windowlength = 5; %you can change that, but keep it odd
hw = (windowlength-1)/2; %hw stands for "half window"
%pre-assigning the output RMS signal
sig_RMS = zeros(1,length(sig));
%I'm adding zeros the the signal's left and right for boundary calculation
sig_zero_padded = [zeros(hw,1), sig, zeros(hw,1)];

sig_zero_padded = sig_zero_padded.^2; %squaring the samples
for i=(hw+1):(length(sig_zero_padded)-hw)
   % calculating square-mean inside the window and assigning the result
   %to the appropriate sample in sig_RMS
  sig_RMS(i-hw) = sqrt(mean(sig_zero_padded(i-hw:i+hw)));
end

Then you set a minimum value for which superior values of the RMS signal will be assigned to 1, and the other to zero :
trigger = 0.5; %I chose arbitrarily
square_sig = (sig_RMS > trigger);

And you should get your square signal.
What we are doing here, in simple language, is taking little bits of the signal and calculating the RMS value of each bit. If this calculated value is significant enough to be considered as "not noise", we assign '1' to it.
EDIT :
I managed to compute it and it worked. Note that I had to take Fs = 10^5 because my computer does not have enough capacity for Fs = 10^6.
Here is the signal (blue), and the calculated RMS signal (red) for windowlength = 8191 :

Then I set trigger = 0.4 and obtained two crystal clear zones identifying the useful signal :

So in conclusion :

Set windowlength to get the RMS signal mean value far enough of the noise while keeping accuracy. In my case it was approximately 1 % of the total number of samples.
Set trigger between the RMS crests and the noise

However, I'm having that feeling that I'm reinventing something so feel free to look up for some literature about this method, I'm sure someone did the math to obtain optimal values of windowslength and trigger. ;)


Matlab: Separate data samples from signal I have data samples for an EMG signal of three arm reps. In the data I have 1 signal that is flat then fluctuates in voltage then becomes flat again shows three pulses for each rep. I want to separate the three reps from the signal (chop the signal into 3 parts). I can manually do this by plotting the signal then zooming in and finding the index before the variance in voltage fluctuates and after it flat lines again. I wanted to know if there is an automated way to do this. Some sort of greedy algorithm or a built in Matlab function. Any ideas?

One example:

[12,-23,1,-2,-40,-45,80,85,14,-28,10,-22,-100,-105,80,85,14,-2,
 16,-2,-200,-45,305,85,1,-2,1,-2]

So I would want:

1 [-40,-45,80,85]
2 [-100,-105,80,85]
3 [-200,-45,305,85]

Usually the flat line values vary from zero.

This is a little hacky, but it uses the three models you have to find an "eigen"-pulse and then takes successive inner products of a normalized signal segment with this "eigen"-pulse to give a measure between -1.0 and 1.0.

The red line in the image is your data.  The black line is this normalized measure.  You'll have to sanity check it (no detections within 2 samples of a previous detection ?), but other than that selecting anything above 0.8 should let you detect the start of the pulses.

R Code Below

  #31201

  arm <- c(12,-23,1,-2,-40,-45,80,85,14,-28,10,-22,-100,-105,80,85,14,-2, 16,-2,-200,-45,305,85,1,-2,1,-2)

  #1 [-40,-45,80,85]
  #2 [-100,-105,80,85]
  #3 [-200,-45,305,85]
  rep1 <- c(-40,-45,80,85)
  rep2 <- c(-100,-105,80,85)
  rep3 <- c(-200,-45,305,85)

  repMat <- rep1 %*% t(rep1) + rep2 %*% t(rep2) + rep3 %*% t(rep3);
  mean1 <- (rep1+rep2+rep2)/3

  eigs <- eigen(repMat)

  plot(rep1, ylim=c(-350, 350), col="red", lwd=6, type="l")
  lines(rep2, col="green", lwd=4)
  lines(rep3, col="blue", lwd=2)
  lines(mean1)
  lines(eigs$vectors[,1]*sqrt(eigs$values[1]), col="purple", lwd=1)

  eigs$vectors[,1] %*% rep1

  similarity <- rep(0,length(arm))
  for (idx in seq(1,length(arm)-3))
  {
    patch <- arm[idx + seq(0,3)]
    patch <- patch / sqrt(sum(patch ^ 2))
    similarity[idx] <- (eigs$vectors[,1] %*% patch)
  }

  plot(similarity, type="l")
  lines(arm/max(arm), col="red")


Path loss effect on BER Does including log distance path loss to a communication system cause a drastic difference to the BER of the system?

In other words, if I have a simulation for a communication system without path loss, does the BER vs SNR curve differ drastically when I add path loss or should both the BER vs SNR curves(with and without path loss) look similar?

Generally speaking, the BER only depends on the SNR. If the received signal power $P_\mathrm{S}$ is reduced due to path loss, while the noise power $P_\mathrm{N}$ remains constant, then the SNR will be reduced, because the SNR is defined as
$$
\gamma = \frac{P_\mathrm{S}}{P_\mathrm{N}}
$$
Consequently, the BER will be increased and the answer to your first question is: yes. However, this should not influence your BER vs SNR curves. Thus the answer to your second question is: no.

If your simulation behaves differently, i.e. if introducing a path loss "shifts" your BER curves, then the SNR is probably not calculated correctly. I suggest that you estimate the received signal power $\hat P_\mathrm{S}$ at the receiver (just before the sampler). Then the power of the additive noise should be calculated by
$$
P_\mathrm{N}=\frac{\hat P_\mathrm{S}}{\gamma}.
$$
Then add noise with mean power $P_\mathrm{N}$ to the received signal. With this setup, your BER/SNR curves should be independent of path loss.


Detecting Central Axis (petiole/rachis) of a Pinnated Palm I'm new in StackExchange and I'd appreciate some help. I'm working in an algorithm for detecting "aguajes", a particular kind of palm. So far, I got good results using a texture feature extractor and other characteristics, but the main problem is that I can't difference between an "aguaje" and other palms (mainly pinnated palms) because the texture is very similar.

1.AGUAJE

2.Pinnated Palm

I see that the main difference in the second one is the presence of a Central Axis (also called Petiole or Rachis). The question is "How to detect it?". I mean, if I'm able to detect the presence of a Petiole, I'd be able to distinguish between this two types of palm. I think color is not an important issue because even when there's no color (gray image), we can discern where is located that Axis.

Thanks in advance.

In your provided image the central axis is straight and have the same color throughout, therefore a Hough operator can detect this.


UART signal recovery from PCM data I have a PCM signal recording of UART data after high-pass filter. The recording looks like this:

I know baud rate and parity setting beforehand. What algorithm could help me recover actual binary data from the recording automatically?

This signal looks pretty good, so that a simple sign decision gives you "low" or "high" periods.
Then, you just start at the beginning, and count samples until you find your start bit.
From the start bit, you just sample your signal in 1 bit period, 2 bit periods, ... and get your transmit data stream. Remember, UART is by definition unsynchronized.


Difference between MATLAB hough and my implementation I am trying to recreate MATLAB's hough function with mine. My code follows

function [H,T,R] = my_hough(x,dr,dtheta)
    rows = size(x,1);
    cols = size(x,2);
    D = sqrt((rows - 1)^2 + (cols - 1)^2);
    Nr = 2*(ceil(D/dr)) + 1;
    diagonal = dr*ceil(D/dr);
    R = -diagonal:dr:diagonal;
    T = -90:dtheta:90-dtheta;
    Ntheta = length(T);
    H = zeros(Nr,Ntheta);
    for i = 1:Ntheta
        for n1 = 1:rows
            for n2 = 1:cols
                if x(n1,n2)==1
                    r = n2*cos(T(i)*pi/180) + n1*sin(T(i)*pi/180);
                    [~,j] = min(abs(R-ones(1,Nr)*r));
                    H(j,i) = H(j,i) + 1;
                end
            end
        end
    end
end

where dr and dtheta are distance and angle resolution. Printing the difference between my Hough table and MATLAB's there are many zeros, but there are also some non-zero elements. Any idea why this is happening?

I don't have the rep to comment on Dima's answer, but the code in the question does convert the $\theta$ to radians already.

This is my implementation of the routine for detecting lines:

[ theta,rho,height,width,Accumulator] = SHTanalyse( m,delta_theta,delta_rho )
%% Get dimensions of input image
[height, width] = size(m);

% Create array for theta
theta = -90:delta_theta:90;
costheta = cosd(theta); % Create lookup tables to optimize computation of cos theta and sin theta
sintheta = sind(theta);

% Find bounds on rho
rhomax = round(sqrt(height^2 + width^2));

% Create array for rho
rho = -rhomax:delta_rho:rhomax;

% Get edgel coordinates
[yEdges,xEdges] =  find(m);

% Initialize the accumulator
Accumulator = zeros(numel(rho),numel(theta));

% Voting process
for k = 1:numel(yEdges) % Loop through edge pixels
    for t = 1:numel(theta)
        % Compute the corresponding rho
        rhvalue = round(xEdges(k).*costheta(t) + yEdges(k).*sintheta(t));
        % Increment accumulator
        Accumulator(rhvalue + round(0.5*numel(rho)) + 1,theta(t) + 90 + 1) = Accumulator(rhvalue + round(0.5*numel(rho)) + 1,theta(t) + 90 + 1) + 1;
    end
end

You can compare its performance with the MATLAB's own Hough function to figure out the discrepancy.


Assess quality of a bandpass filter without human intervention I wanted to know if there is a way to assess the quality of a bandpass filter created by the Remez exchange algorithm. I am using the Remez algorithm to determine the best bandpass filter at a number of different sample frequencies and numbers of coefficients and was wondering if there some easy way to assess the quality of the bandpass filter created without just looking at it's frequency response?

I want to be able to iteratively explore the parameter space and return the sample frequency and number of coefficients that produce the best band-pass filter. I am aware that there is not an objective "best" filter, by "best" I mean the filter that best matches a certain specification.

The Remez exchange algorithm not only returns the optimal filter coefficients but also the approximation error $\delta$. This error determines the maximum amplitude error in the pass bands and in the stop bands. The actual amplitude error in band $i$ is given by $\delta/W_i$, where $W_i$ is the error weight in band $i$.

If the algorithm converges, then you know for sure that the maximum errors $\delta/W_i$ occur at the band edges of band $i$. And you know the band edges because they are part of the specs.


Calculate bag of words feature vector In visual bag of words model, I have been able to construct the visual codebook through kmeans clustering of SIFT descriptors. How to calculate the feature vector for an image then?

P/S:

For each image, we can find interesting SIFT points, and for each points we have a SIFT descriptor (which is usually a 128 length vector).

im1 ==> SIFT feature f1 (10 by 128) (here 10 is an abitrary number)
im2 ==> SIFT feature f2 (20 by 128)
...

If we combine all SIFT features, f=[f1; f2; ..] and perform kmeans clustering we will get the codebook c=[c1; c2; .. c10] which is bow codebooks.

From the codebook how can we find the feature vectore, represent image im1?

The feature vector is just the histogram of how many times a feature from each cluster appeared in the image.


How can I shift this signal to it's midpoint in order to decode it I'm wondering if there's a simple way to decode this pulse signal as I'm reading it in. I'm thinking the simplest way to decode this is if i can shift the whole signal down so that it's centered around 0. Then if my pulse goes above 0, I know it's a 1 and it'f it's below zero, i know it's a zero. Figures below to see what i'm trying to do. Please note that my sampling frequnecy $Fs = 8Hz$, my pulse width is between 1-2 seconds, and come in pairs of 2. these pairs are spaced about 20 seconds apart. There is also a 3 pulse header which denotes the start of the message, with the first pulse being wider than the other pulses.

if I know the peak value of the pulses, then this is obviously just taking my signal and subtracting half of the peak, $d(t) = r(t) - peak/2$. The problem is this peak amplitude is going to vary, as is the noise level. If i try to just subtract a running average, then the noise part of my signal will be centered around 0 as well, which then my decision won't be a simple above/below zero decision. I've plotted subtracting the mean in figure 3 below. My signal is $d(t) = r(t) - runAvg(past 4 samples)$

It seems like there should be an easy way to filter this, either in the time or frequency domain, I'm just not thinking of it. Thanks for your time!

You don't need to remove the DC part of your signal (i.e. "to center" it). In terms of decoding performance, there is no gain. Your signal looks like a typical ASK. Applying a simple integrator and sampling at $t=T$, where $T$ is the pulse duration, would do the job.

You only need to define a threshold, which usually corresponds to $AT/2$, where $A$ is the expected pulse amplitude.

Edit: forgot to mention the fact that you need to determine when to sample the output of the integrator. A symbol-level synchroniser is required for that.


(How to ask a Homework Question): Define poles by using proportional controller Given is a process with the transfer function

$$G(s) = \frac{s - 1}{s^2 + 3s + 2}$$

I want to create a controller so that the poles of the controlled system are

$$p_{1,2} = -4 \pm i$$

Is it possible to do this using a proportional controller?

I already calculated the transfer function of the closed loop system, which is:

$$G_1(s) = \frac{k \cdot (s-1)}{(s+1) \cdot (s+2) + k \cdot (s-1)}$$

Then I tried to figure out whether there exists a value for k which leads to the wanted poles and saw that there exists no such value, so my answer would be that it is not possible to use a proportional controller. Please correct me if my solution is wrong.

Nice work and your conclusion is correct.

Further confirmation is the poles in the closed loop system are the roots of $$1+K G(s)=0$$

Which results in the polynomial with a numerator

$$s^2+(3+K)s+(2-K) $$

In order to have complex roots, using the quadratic formula $\sqrt{b^2-4ac}$ must be negative.

So we want to test if $4(2-K)>(3+K)^2$ for any K

Since $4(2-K)=(3+K)^2$ has roots at approximately -.101, -9.898:

$$K^2+10K+1=0$$

The range of K where the poles can be complex is therefore from -.101 to -9.898.

Showing all poles for the closed loop system as a function of K (root locus) over the range of K=-10.5 to 0 as plotted below confirms that $-4 \pm i$ is not a closed loop pole.

We could mathematically prove this by using the relationship of "Break-Outs" on the root locus:

(See Dr. Cheever's excellent pages at http://lpsa.swarthmore.edu/Root_Locus/Example2/Example2.html)

As Dr. Cheever summarized, break-out or break-in points occur where N(s)D'(s)-N'(s)D(s)=0 where N(s) and D(s) are the numerator and denominator polynomial of the open loop transfer function, and N'(s) and D'(s) are the differentiation.

$N(s)= s-1$

$N'(s) = 1$

$D(s) = s^2+s3+2$

$D'(s) = 2s+3$

Solving this results in roots at approximately -1.4495 and 3.4495 indicating the break out and break in locations. The root locus will break out at -1.4495 when K=-0.101 and follow a trajectory as shown in the diagram above to the break-in at 3.4495 when K=-9.898. Therefore $-4 \pm i$ is not a closed loop pole.

Note: Root locus was plotted using Octave with the rlocus command as part of the control toolbox

pkg load control %if not already loaded
sys=tf([1 -1],[1 3 2]);
rlocus(sys,.1,-10.5,0);


How to perform convolution using polyphase structure I am trying to implement convolution using polyphase structure but, how to split input samples among different filter banks is my major doubt. Please explain the same with example?

I am not quite sure, if i understand you question properly, it would be better if you ask with detail and maybe with a figure, nevertheless if you want convolution of co-efficents of the transfer function then below might be of help. :

and with regard to splitting.. you can easily do it with something like this :

x(1:2:end) = cos(2*pi*f*t(1:2:end));
x(2:2:end) = cos(2*pi*f*t(2:2:end));


Confusion regarding sinusoidal synthesis formula, I am reading signal processing first by McClellan

As shown in attached snapshot, we have formula of sinusoidal synthesis

I am confused ,why in expression of x(t), there is a '2' in first line immediately after summation sign but that'2' disappear in 2nd line?

Simple:

Euler's formula:

$$ e^{i\theta} = \cos(\theta) + i \sin(\theta) $$

and

$$ e^{-i\theta} = \cos(\theta) - i \sin(\theta) $$

Add them together:

$$ e^{i\theta} + e^{-i\theta} = 2 \cos(\theta) $$

You should see that in your equation.

From there it is usually stated:

$$ \cos(\theta) = \frac{ e^{i\theta} + e^{-i\theta}}{ 2 } $$

See my article

The Exponential Nature of the Complex Unit Circle

for a more expansive explanation.


There appears to be a contradiction with the recurrence solutiom Consider the equation :
$$
y[n]-\frac{1}{2}y[n-1]=x[n]
$$
where $x[n]:=\left(\frac{1}{3}\right)^{n}u[n]$. First, I founded the homogenous solution $y_{h}:=A\frac{1}{2^{n}}$
$$
y_{h}[n]=\frac{1}{2^{n}}y[0]
$$
where I insist $A=y[0]$. Furthermore, The particular solution is given of the form $y_{p}[n]:=B\left(\frac{1}{3}\right)^{n}$ I shall determine $B$ as follow :
$$
B\left(\frac{1}{3}\right)^{n}u[n]-\frac{1}{2}B\left(\frac{1}{3}\right)^{n-1}u[n-1]=\left(\frac{1}{3}\right)^{n}u[n]
$$
\begin{align*}
\implies B&=\frac{\displaystyle\left(\frac{1}{3}\right)^{n}u[n]}{\displaystyle\left(\frac{1}{3}\right)^{n}u[n]-\frac{1}{2}\displaystyle\left(\frac{1}{3}\right)^{n-1}u[n-1]}\\ \\
&=\frac{\displaystyle\left(\frac{1}{3}\right)^{n}}{\displaystyle\left(\frac{1}{3}\right)^{n}-\displaystyle\frac{1}{2}\left(\frac{1}{3}\right)^{n-1}}\\ \\
&=-2
\end{align*}
So now that we have $B=-2$, We are given the initial rest condition so what I did was sub $n=0$ in difference equation
$$
y[0]-\frac{1}{2}y[-1]=x[0]
$$
$$
\implies y[0]=1
$$
$$
**\text{So this means $A=y[0]=1$}**
$$
However, the major issue is that :
$$
y[0]-\frac{1}{2}y[-1]=x[0]
$$
$$
\implies y[0]=1
$$
$$
\implies y[n]:= y_{h}[n]+y_{p}[n] \implies y[0]=y_{h}[0]+y_{p}[0]=1
$$
$$
\implies A\left(\frac{1}{2}\right)^{0} -2\left(\frac{1}{3}\right)^{0}=1
$$
$$
\implies A=3
$$
So I have $A=1$ and $A=3$, an obvious contradiction and I hope someone can help me

From the comments under Tim Wescott's answer, it seems that there is still  some confusion concerning the correct way to solve such problems. Let's consider the non-homogeneous difference equation (DE)
$$y[n]-p y[n-1]=q^n,\qquad p\neq q\tag{1}$$
with initial value $y[-1]=c$. Note that initially we don't restrict the forcing function on the right-hand side of $(1)$ to equal zero for $n<0$.
As correctly shown in the OP, the general solution to the homogeneous DE has the form
$$y_h[n]=Ap^n\tag{2}$$
which is also called the natural response,
and a particular solution (forced response) of $(1)$ is given by
$$y_p[n]=Bq^n\tag{3}$$
The constant $B$ is determined by substituting $(3)$ into $(1)$:
$$Bq^n-pBq^{n-1}=q^n\tag{4}$$
from which we obtain $B=\frac{q}{q-p}$. Consequently, the general solution of the non-homogeneous DE $(1)$ is
$$y[n]=y_h[n]+y_p[n]=Ap^n+\frac{q}{q-p}q^n\tag{5}$$
Note that the solution $(5)$ is valid for all values of $n$.
The value of $A$ is determined by the initial condition $y[-1]=c$. Substituting $n=-1$ into $(5)$ and setting it equal to $c$ results in
$$A=cp-\frac{p}{q-p}\tag{6}$$
The solution to the initial value problem is hence given by
$$y[n]=cp^{n+1}+\frac{1}{q-p}\left[q^{n+1}-p^{n+1}\right]\tag{7}$$
We now only consider values $n\ge 0$, simply because that's the range we're interested in.
Note that the term $cp^{n+1}$ in $(7)$ is the zero-input response (ZIR), and the remaining term is the zero-state response (ZSR).
From an engineering point of view, composing the solution into a ZSR and a ZIR appears more natural and more useful than the classical decomposition into a homogeneous and a particular solution, because the decomposition of the solution into ZSR and ZIR separates the contributions arising from an external input and from internal conditions, respectively.


Shouldn't the Sampling Theorem imply that there should be no information loss at all after a signal is processed? I am very, very new to signal processing (only started a few days ago) so please bear with me because I may be missing the big picture.

Consider an arbitrary signal $x(t)$ with varying frequency to be processed by a computer.

Given the right sampling rate w.r.t the frequency, shouldn't the Sampling Theorem essentially imply that no information will be lost at all after the sampled signal is converted back to its continuous form, since there will be a sufficient number of discrete points for which the signal is described?

For no information to be lost on conversion back to continuous form, the signal would first need to be perfectly band-limited, and you would need an ideal reconstruction filter.  A perfectly band-limited signal is infinite in extent. Since you want this arbitrary signal to be processed by a computer, your computer would need infinite memory.  You would also have to wait an infinite amount of time for your perfect reconstruction filter to settle.

Finite sampling time and filter length problems also completely ignores other potential lossy issues, such as quantization and sampling clock jitter, etc.

One might be able to process an infinite signal that can be represented in certain closed symbolic forms within finite time and resources, but there does not seem to be a method to convert any arbitrary signal into such a form.

Back in the real-world (stuff you can buy or make), one would normally accept an information loss from imperfect band-limiting, finite filter length, jitter (etc.) that is around or below the quantization and numeric noise floors.  This allows the processing to happen using an amount of RAM that one can afford to hopefully finish within one's lifetime.  Thus leading to information loss and imperfect reconstruction.